C///////////////////////////////////////////////////////////////////////
C
C     TWOPNT
C
C     TWO POINT BOUNDARY VALUE PROBLEM SOLVER
C
C     THIS IS VERSION 2.22B OF MARCH 1991.
C
C     WRITTEN BY: DR. JOSEPH F. GRCAR
C                 SANDIA NATIONAL LABORATORY
C                 SCIENTIFIC COMPUTING DIVISION 8233
C                 LIVERMORE, CALIFORNIA  94551-0969  USA
C
C                 (415) 294-2662
C                 (FTS) 234-2662
C
C///////////////////////////////////////////////////////////////////////
C
C     CHANGES FROM THE PREVIOUS VERSION:
C
C     1) ALLOW NEGATIVE STEPS0 AS FLAG TO AVOID NEWTON SEARCH.
C
C///////////////////////////////////////////////////////////////////////
C
C     ARGUMENTS FOR SUBROUTINE TWOPNT:
C
C     ERROR   OUTPUT LOGICAL - ERROR FLAG.  IF TRUE, THEN A PROGRAMMING
C             ERROR OR A SEVERE NUMERICAL ERROR BLOCKS EXECUTION.  ERROR
C             MESSAGES APPEAR IN THE OUTPUT TEXT FILE (NO MATTER THE
C             OUTPUT LEVEL).
C
C     TEXT    INPUT INTEGER - UNIT NUMBER FOR AN OUTPUT FILE.  ALL TEXT
C             OUTPUT GOES TO THIS UNIT.  'LEVEL' CONTROLS THE QUANTITY
C             OF OUTPUT, BUT ZERO AND NEGATIVE VALUES FOR 'TEXT' SUSPEND
C             ALL OUTPUT (INCLUDING ERROR MESSAGES).
C
C     LEVEL   INPUT INTEGER DIMENSIONED (2) - OUTPUT LEVEL.  IN THE TREE
C             OF SUBROUTINE CALLS STEMMING FROM TWOPNT, 'LEVEL(1)'
C             SPECIFIES THE HIGHEST ROUTINES THAT WRITE INFORMATIVE
C             MESSAGES TO THE OUTPUT TEXT FILE.  IF 'LEVEL(1)' EQUALS 1,
C             THEN ONLY TWOPNT WRITES MESSAGES; IF 'LEVEL(1)' EQUALS 2,
C             THEN SUBROUTINES CALLED DIRECTLY BY TWOPNT ALSO WRITE
C             MESSAGES; AND SO ON.  IF 'LEVEL(1)' IS ZERO OR NEGATIVE,
C             THEN ONLY ERROR MESSAGES APPEAR.  IF THE UNIT NUMBER
C             'TEXT' IS ZERO OR NEGATIVE, THEN NOT EVEN ERROR MESSAGES
C             APPEAR.
C
C             'LEVEL(2)' SPECIFIES THE HIGHEST ROUTINES THAT WRITE THE
C             SOLUTION TO THE OUTPUT FILE.  IF 'LEVEL(2)' EQUALS 1, THEN
C             THE SOLUTION APPEARS AT THE START AND FINISH OF TWOPNT; IF
C             'LEVEL(2)' EQUALS 2, THEN THE SOLUTION APPEARS AT THE
C             START OF TWOPNT AND AT THE FINISH OF EVERY SUBROUTINE
C             CALLED DIRECTLY BY TWOPNT (THE SOLUTION DOES NOT REAPPEAR
C             WHEN TWOPNT ITSELF FINISHES); AND SO ON.  'LEVEL(2)' MUST
C             BE LESS THAN OR EQUAL TO 'LEVEL(2)'.
C
C     ISIZE   INPUT INTEGER - SIZE OF THE INTEGER WORK SPACE ARRAY
C             'IWORK'.  MUST BE AT LEAST 3 * 'PMAX'.
C
C     IWORK   INTEGER DIMENSIONED (ISIZE) - INTEGER WORK SPACE.
C
C     RSIZE   REAL INTEGER - SIZE OF THE REAL WORK SPACE ARRAY 'RWORK'.
C             MUST BE AT LEAST 7 * 'SCALRS' + (7 * 'COMPS' + 2) *
C             'PMAX'.
C
C     RWORK   REAL DIMENSIONED (ISIZE) - REAL WORK SPACE.
C
C     ABOVE   INPUT/OUTPUT REAL DIMENSIONED (SCALRS + COMPS * PMAX) -
C             UPPER BOUNDS.  EACH SOLUTION VARIABLE HAS ITS OWN BOUND.
C             NO SOLUTION VARIABLE EXCEEDS ITS UPPER BOUND AT ANY TIME.
C             IF MESH REFINEMENT OCCURS, THEN NEW POINTS RECEIVE BOUNDS
C             INTERPOLATED FROM SURROUNDING POINTS.  THE USER MAY
C             OVERWRITE THESE DURING 'UPDATE' REVERSE COMMUNICATION
C             REQUESTS.
C
C     ACTIVE  INPUT LOGICAL DIMENSIONED (COMPS) - MARKERS FOR ACTIVE
C             COMPONENTS.  ONLY ACTIVE COMPONENTS OF THE SOLUTION GOVERN
C             MESH REFINEMENT.  COMPONENT J IS ACTIVE WHEN 'ACTIVE(J)'
C             IS TRUE.
C
C     ADAPT   INPUT LOGICAL - REFINEMENT DIRECTIVE.  IF TRUE, THEN
C             PERMIT MESH REFINEMENT.
C
C     BELOW   INPUT/OUTPUT REAL DIMENSIONED (SCALRS + COMPS * PMAX) -
C             LOWER BOUNDS.  EACH SOLUTION VARIABLE HAS ITS OWN BOUND.
C             NO SOLUTION SOLUTION VARIABLE FALLS SHORT OF ITS LOWER
C             BOUND AT ANY TIME.  IF MESH REFINEMENT OCCURS, THEN NEW
C             POINTS RECEIVE BOUNDS INTERPOLATED FROM SURROUNDING
C             POINTS.  THE USER MAY OVERWRITE THESE DURING 'UPDATE'
C             REVERSE COMMUNICATION REQUESTS.
C
C     BINARY  INPUT INTEGER - UNIT NUMBER FOR AN OUTPUT FILE.  THIS UNIT
C             RECEIVES BINARY OUTPUT THAT SUMMARIZES THE PROGRESS OF
C             TWOPNT.  PROGRAM 'SHOW' GRAPHS THE DATA.  ZERO AND
C             NEGATIVE VALUES FOR 'BINARY' SUSPEND ALL OUTPUT.
C
C     BUFFER  INPUT/OUTPUT REAL DIMENSIONED (SCALRS + COMPS * PMAX) -
C             REVERSE COMMUNICATION BUFFER.
C
C     COMPS   INPUT INTEGER - NUMBER OF COMPONENTS.
C
C     CONDIT  INPUT REAL - CONDITION NUMBER.  'CONDIT' PASSES THE
C             CONDITION NUMBER OF THE LATEST JACOBIAN MATRIX.  ZERO AND
C             NEGATIVE VALUES MEAN THE NUMBER IS UNKNOWN.
C
C     FUNCT   OUTPUT LOGICAL - REVERSE COMMUNICATION SIGNAL.  IF TRUE,
C             THEN (1) EVALUATE THE RESIDUAL FUNCTION AT THE APPROXIMATE
C             SOLUTION IN 'BUFFER', (2) STORE THE RESULT IN 'BUFFER',
C             AND (3) REENTER TWOPNT.  WHEN 'TIME' IS ALSO TRUE, THEN
C             THE RESIDUAL FUNCTION SHOULD BE THE ONE FOR THE
C             TIME-DEPENDENT PROBLEM.
C
C     JACOB   OUTPUT LOGICAL - REVERSE COMMUNICATION SIGNAL.  IF TRUE,
C             THEN (1) EVALUATE THE JACOBIAN MATRIX OF THE RESIDUAL
C             FUNCTION AT THE APPROXIMATE SOLUTION IN 'BUFFER', (2)
C             PREPARE TO SOLVE SYSTEMS OF LINEAR EQUATIONS HAVING THIS
C             MATRIX OF COEFFICIENTS (FOR EXAMPLE, FACTOR THE MATRIX),
C             (3) OPTIONALLY PLACE THE CONDITION NUMBER OF THE MATRIX IN
C             'CONDIT', AND (4) CALL TWOPNT AGAIN.  WHEN 'TIME' IS ALSO
C             TRUE, THEN THE RESIDUAL FUNCTION SHOULD BE THE ONE FOR THE
C             TIME-DEPENDENT PROBLEM.
C
C     MARK    OUTPUT LOGICAL DIMENSIONED (PMAX) - MARKERS FOR NEW
C             POINTS.  IF 'UPDATE' IS TRUE, THEN 'MESH(K)' IS NEW
C             EXACTLY WHEN 'MARK(K)' IS TRUE.
C
C     MESH    INPUT/OUTPUT REAL DIMENSIONED (PMAX) - MESH POINTS.  THIS
C             ARRAY CONTAINS THE VALUES OF THE INDEPENDENT SPATIAL
C             VARIABLE AT WHICH THE STEADY-STATE AND TIME-DEPENDENT
C             FUNCTIONS DISCRETIZE THE DIFFERENTIAL EQUATIONS.  THE USER
C             CHOOSES THE INITIAL POINTS.  THE POINTS MUST LIE IN
C             INCREASING OR DECREASING ORDER.  IF 'ADAPT' IS TRUE, THEN
C             TWOPNT MAY ADD POINTS UP TO THE CAPACITY OF THE ARRAY.
C             'POINTS' IS THE CURRENT NUMBER.
C
C     PASS    OUTPUT INTEGER - NUMBER OF THE LATEST PROBLEM PASS ON THE
C             LATEST MESH.  SEE 'PASSES'.
C
C     PASSES  INPUT INTEGER - NUMBER OF PASSES PER MESH.  IF IT IS
C             DESIRABLE TO POSE SEVERAL PROBLEMS (PERHAPS OF INCREASING
C             DIFFICULTY) ON EACH MESH, THEN SET PASSES TO THE NUMBER OF
C             PROBLEMS.  OTHERWISE, SET PASSES TO ONE.
C
C     PLIMIT  INPUT INTEGER - MAXIMUM NUMBER OF POINTS THAT MAY BE
C             SIMULTANEOUSLY ADDED TO THE MESH.
C
C     PMAX    INPUT INTEGER - MAXIMUM NUMBER OF MESH POINTS.  PMAX
C             DETERMINES THE SIZE OF MANY ARRAYS.
C
C     POINTS  INPUT/OUTPUT INTEGER - CURRENT NUMBER OF MESH POINTS.  THE
C             COMPONENTS FOR EACH POINT LIE TOGETHER AFTER THE SCALAR
C             VARIABLES IN ARRAYS 'ABOVE', 'BELOW', 'BUFFER', AND 'X'.
C
C     REENT   INPUT/OUTPUT LOGICAL - REENTER SIGNAL.  SET FALSE ON INPUT
C             TO MARK THE START OF A NEW PROBLEM.  SET TRUE ON OUTPUT
C             WHEN TWOPNT REQUESTS INFORMATION.  IN THIS CASE, PERFORM
C             THE TASK INDICATED BY ONE OF THE REVERSE COMMUNICATION
C             SIGNALS ('FUNCT', 'JACOB', 'SAVE', 'SHOW', 'SOLVE',
C             'STORE', OR 'UPDATE') AND REENTER TWOPNT.
C
C     REPORT  OUTPUT INTEGER - FAILURE NOTE.  IF BOTH 'ERROR' AND
C             'SUCCES' ARE FALSE, THEN 'REPORT' EXPLAINS WHY (OTHERWISE
C             'REPORT' RETURNS ZERO).  (1) NO CHANGE:  THE NEWTON
C             ALGORITHM DID NOT SOLVE THE STEADY-STATE PROBLEM, AND THE
C             SUBSEQUENT TIME STEP DID NOT CHANGE THE SOLUTION.  THE
C             TIME STEP IS PROBABLY TOO SMALL.  (2) NO SOLUTION:  THE
C             NEWTON ALGORITHM DID NOT SOLVE A TIME-DEPENDENT PROBLEM.
C             THE TIME STEP IS PROBABLY TOO LARGE.  (3) NO SPACE:  THE
C             MESH REFINEMENT CRITERIA ARE NOT SATISFIED ON THE LARGEST
C             POSSIBLE GRID.
C
C     SAVE    OUTPUT LOGICAL - REVERSE COMMUNICATION SIGNAL.  IF TRUE,
C             THEN THE LATEST APPROXIMATE SOLUTION MAY BE USED AS A
C             STARTING ESTIMATE SHOULD TWOPNT HAVE TO BE RESTARTED (DUE
C             TO SYSTEM CRASH OR TIME OUT, FOR EXAMPLE).  THUS, (1)
C             OPTIONALLY SAVE 'MESH', 'POINTS', AND THE SOLUTION FOUND
C             IN 'BUFFER', AND (2) REENTER TWOPNT.
C
C     SCALRS  INPUT INTEGER - NUMBER OF SCALAR VARIABLES.  THESE PRECEDE
C             THE MESH VARIABLES IN THE ARRAYS X AND BUFFER.
C
C     SHOW    OUTPUT LOGICAL - REVERSE COMMUNICATION SIGNAL.  IF TRUE,
C             THEN (1) WRITE TO UNIT 'TEXT' THE APPROXIMATE SOLUTION IN
C             'BUFFER', AND (2) REENTER TWOPNT.  THE USER CAN DISPLAY
C             THE SOLUTION WITH PROBLEM DEPENDENT DATA, SUCH AS
C             COMPONENT LABELS, THAT ARE UNKNOWN TO TWOPNT.  'LEVEL(2)'
C             DETERMINES WHEN TWOPNT REQUESTS THE SOLUTION BE SHOWN.
C
C     SOLVE   OUTPUT LOGICAL - REVERSE COMMUNICATION SIGNAL.  IF TRUE,
C             THEN (1) SOLVE A SYSTEM OF LINEAR EQUATIONS HAVING THE
C             RIGHT HAND SIDE IN 'BUFFER' AND USING THE MOST RECENT
C             JACOBIAN MATRIX, (2) STORE THE SOLUTION IN 'BUFFER', AND
C             (3) REENTER TWOPNT.
C
C     SSABS   INPUT REAL - ABSOLUTE BOUND FOR THE STEADY-STATE PROBLEM.
C             THE NEWTON ALGORITHM TERMINATES IF EACH ENTRY OF THE NEXT
C             (UNDAMPED) STEP IS EITHER ABSOLUTELY SMALL OR SMALL
C             RELATIVE TO THE CORRESPONDING ENTRY OF THE APPROXIMATE
C             SOLUTION.
C
C     SSAGE   INPUT INTEGER - AGE OF THE STEADY STATE JACOBIAN.  NUMBER
C             OF NEWTON STEPS TO BE TAKEN BEFORE RECOMPUTING THE
C             JACOBIAN.
C
C     SSREL   INPUT REAL - RELATIVE BOUND FOR THE STEADY-STATE PROBLEM.
C             THE NEWTON TERMINATES IF EACH ENTRY OF THE NEXT (UNDAMPED)
C             STEP IS EITHER ABSOLUTELY SMALL OR SMALL RELATIVE TO THE
C             CORRESPONDING ENTRY OF THE APPROXIMATE SOLUTION.
C
C     STEPS0  INPUT INTEGER - NUMBER OF INITIAL TIME STEPS.  IF
C             POSITIVE, THEN THIS MANY STEPS SHOULD BE COMPLETED OF THE
C             TIME-DEPENDENT PROBLEM BEFORE ATTEMPTING THE NEWTON
C             ALGORITHM.  IF NEGATIVE, NEWTON IS NOT ATTEMPTED.
C
C     STEPS1  INPUT INTEGER - NUMBER OF SMOOTHING TIME STEPS.  IF THE
C             NEWTON ALGORITHM CANNOT SOLVE THE STEADY-STATE PROBLEM,
C             THEN AT MOST THIS MANY STEPS SHOULD BE ATTEMPTED OF THE
C             TIME-DEPENDENT PROBLEM.
C
C     STEPS2  INPUT INTEGER - NUMBER OF TIME STEPS TO BE TAKEN BEFORE
C             INCREASING THE TIME STEP.
C
C     STORE   OUTPUT LOGICAL - REVERSE COMMUNICATION SIGNAL.  IF TRUE,
C             THEN TWOPNT IS ABOUT TO PERFORM OR IS PERFORMING A TIME
C             INTEGRATION.  'BUFFER' CONTAINS EITHER THE INITIAL VALUE
C             OR THE SOLUTION FROM THE LAST TIME STEP.  THE
C             TIME-DEPENDENT FUNCTION NEEDS THIS TO EVALUATE TIME
C             DERIVATIVES.  THUS, (1) STORE THE SOLUTION, AND (2)
C             REENTER TWOPNT.
C
C     SUCCES  OUTPUT LOGICAL - COMPLETION STATUS.  IF TRUE, THEN TWOPNT
C             SOLVES THE BOUNDARY VALUE PROBLEM TO THE EXTENT THAT THE
C             SOLUTION IN 'X' SATISFIES THE MESH REFINEMENT CRITERIA.
C             IF FALSE, THEN THE SOLUTION IN 'X' DOES NOT SATISFY THE
C             CRITERIA ('REPORT' MAY EXPLAIN WHY).
C
C     TDABS   INPUT REAL - ABSOLUTE BOUND FOR THE TIME-DEPENDENT
C             PROBLEM.  THE NEWTON ALGORITHM TERMINATES IF EACH ENTRY OF
C             THE NEXT (UNDAMPED) STEP IS EITHER ABSOLUTELY SMALL OR
C             SMALL RELATIVE TO THE CORRESPONDING ENTRY OF THE
C             APPROXIMATE SOLUTION.
C
C     TDAGE   INPUT INTEGER - AGE OF THE TIME DEPENDENT JACOBIAN.
C             NUMBER OF NEWTON STEPS TO BE TAKEN BEFORE RECOMPUTING THE
C             JACOBIAN DURING TIME STEPPING.
C
C     TDEC    INPUT REAL - FACTOR BY WHICH TO DIVIDE THE TIME STEP WHEN
C             NECESSARY.
C
C     TDREL   INPUT REAL - RELATIVE BOUND FOR THE TIME-DEPENDENT
C             PROBLEM.  THE NEWTON ALGORITHM TERMINATES IF EACH ENTRY OF
C             THE NEXT (UNDAMPED) STEP IS EITHER ABSOLUTELY SMALL OR
C             SMALL RELATIVE TO THE CORRESPONDING ENTRY OF THE
C             APPROXIMATE SOLUTION.
C
C     TIME    OUTPUT LOGICAL - REVERSE COMMUNICATION SIGNAL.  IF TRUE,
C             THEN SATISFY THE ACCOMPANYING REQUEST FOR 'FUNCT' OR
C             'JACOB' USING THE TIME-DEPENDENT FUNCTION.
C
C     TINC    INPUT REAL - FACTOR BY WHICH TO MULTIPLY THE TIME STEP
C             WHEN THE NUMBER OF STEPS SPECIFIED BY 'RETIRE' HAS PASSED.
C
C     TMAX    INPUT REAL - MAXIMUM TIME STEP.
C
C     TMIN    INPUT REAL - MINIMUM TIME STEP.
C
C     TOLER0  INPUT REAL - ABSOLUTE AND RELATIVE SIGNIFICANCE LEVEL.
C             THE MESH REFINEMENT CRITERIA ('TOLER1' AND 'TOLER2') APPLY
C             ONLY TO ACTIVE, SIGNIFICANT COMPONENTS.  A COMPONENT IS
C             INSIGNIFICANT WHEN ITS CHANGE (MAXIMUM LESS MINIMUM) IS
C             EITHER ABSOLUTELY SMALL OR SMALL RELATIVE TO ITS MAXIMUM
C             MAGNITUDE.
C
C     TOLER1  INPUT REAL - BOUND UPON THE RANGE OF A COMPONENT.  THE
C             FIRST MESH REFINEMENT CRITERION REQUIRES THAT EACH
C             COMPONENT'S CHANGE OVER ADJACENT MESH POINTS BE NO MORE
C             THAN THE FRACTION 'TOLER1' OF THE COMPONENT'S CHANGE
C             (MAXIMUM LESS MINIMUM) OVER ALL MESH POINTS.  IF SOME
C             COMPONENT VIOLATES THIS CRITERION OVER SOME MESH INTERVAL,
C             THEN TWOPNT ADDS THE INTERVAL'S MIDPOINT TO THE MESH.
C
C     TOLER2  INPUT REAL - BOUND UPON THE RANGE OF A COMPONENT'S
C             DERIVATIVE.  THE SECOND MESH REFINEMENT CRITERION REQUIRES
C             THAT THE CHANGE IN EACH COMPONENT'S DERIVATIVE OVER
C             ADJACENT MESH INTERVALS BE NO MORE THAN THE FRACTION
C             'TOLER2' OF THE DERIVATIVE'S CHANGE (MAXIMUM LESS MINIMUM)
C             OVER ALL MESH INTERVALS.  IF SOME COMPONENT VIOLATES THIS
C             CRITERION OVER SOME PAIR OF ADJACENT MESH INTERVALS, THEN
C             TWOPNT ADDS THE INTERVALS' MIDPOINTS TO THE MESH.
C
C     TSTEP0  INPUT REAL - INITIAL TIME STEP.
C
C     TSTEP1  OUTPUT REAL - LATEST SIZE OF THE TIME STEP TO BE USED IN
C             EVALUATING THE FUNCTION (WHEN TIME STEPS ARE BEING TAKEN).
C
C     UPDATE  OUTPUT LOGICAL - REVERSE COMMUNICATION SIGNAL.  IF TRUE,
C             THEN 'MESH' CONTAINS SOME NEW POINTS (MARKED BY 'MARK').
C             TWOPNT CHOOSES INTERPOLATED BOUNDS AND INTERPOLATED
C             SOLUTION ESTIMATES FOR THE NEW POINTS.  THESE MAY BE
C             INADEQUATE.  THUS (1) OPTIONALLY ADJUST 'ABOVE', 'BELOW',
C             AND THE APPROXIMATE SOLUTION IN 'BUFFER', AND (2) REENTER
C             TWOPNT.
C
C     X       INPUT/OUTPUT REAL DIMENSIONED (SCALRS + COMPS * PMAX) -
C             SOLUTION.  ON INPUT, AN ESTIMATE.  ON OUTPUT, THE
C             SOLUTION.
C
C///////////////////////////////////////////////////////////////////////

      SUBROUTINE COPY (N, X, Y)

      IMPLICIT COMPLEX (A - Z)
C*****precision > double
      DOUBLE PRECISION X, Y
C*****END precision > double
      INTEGER J, N
C*****precision > single
C      REAL X, Y
C*****END precision > single

      DIMENSION X(N), Y(N)

      DO 0100 J = 1, N
         Y(J) = X(J)
0100  CONTINUE

      RETURN
      END
      SUBROUTINE CPUTIM (TIMER)

C///////////////////////////////////////////////////////////////////////
C
C     CPUTIM
C
C     OBTAIN ELAPSED CPU JOB TIME IN SECONDS.
C
C///////////////////////////////////////////////////////////////////////

      IMPLICIT COMPLEX (A - Z)
C*****precision > double
      DOUBLE PRECISION TIMER
C*****END precision > double
C*****precision > single
C      REAL TIMER
C*****END precision > single

C     THE SYSTEM LIBRARY BASELIB CONTAINS SUBROUTINE IZM00.
C*****CRAY/CTSS (LIVERMORE ENVIRONMENT)
C      EXTERNAL IZM00
C      INTEGER FLAG, IZM00, VALUE
C      REAL TEMP
C      DIMENSION VALUE(5)
C
C      FLAG = IZM00 (VALUE)
C      IF (FLAG .EQ. 0) THEN
C         TEMP = REAL (VALUE(1)) / 1.0E6
C      ELSE
C         TEMP = 0.0
C      END IF
C*****END CRAY/CTSS (LIVERMORE ENVIRONMENT)

C     THE SYSTEM LIBRARY CFTLIB CONTAINS SUBROUTINE ISECOND.
C*****CRAY/CTSS (LOS ALAMOS ENVIRONMENT)
C      EXTERNAL ISECOND
C      INTEGER VALUE
C      REAL TEMP
C
C      CALL ISECOND (VALUE)
C      TEMP = REAL (VALUE) / 1.0E6
C*****END CRAY/CTSS (LOS ALAMOS ENVIRONMENT)

C*****CRAY/UNICOS
C      REAL SECOND, TEMP
C      TEMP = SECOND ()
C*****END CRAY/UNICOS

C*****SUN/UNIX
      EXTERNAL REAL ETIME
      REAL LIST
      DIMENSION LIST(2)

      TEMP = ETIME (LIST)
C*****END SUN/UNIX

C     THE VAX/VMS VERSION USES SYSTEM SERVICE CALLS.
C*****VAX/VMS
C      INTEGER*2 LIST
C      INTEGER*4 ADDRESS, TICS
C      REAL*4 TEMP
C
C      DIMENSION LIST(6)
C
C      EQUIVALENCE (LIST(3), ADDRESS)
C
C      DATA LIST(1), LIST(2), LIST(5), LIST(6) / 4, '0407'X, 0, 0 /
C
C      ADDRESS = %LOC (TICS)
C      CALL SYS$GETJPIW (, , , LIST, , , )
C      TEMP = REAL (TICS) / 1.0E2
C*****END VAX/VMS

C*****precision > double
      TIMER = DBLE (TEMP)
C*****END precision > double
C*****precision > single
C      TIMER = TEMP
C*****END precision > single

      RETURN
      END
      SUBROUTINE ELAPSE (TIMER)

      IMPLICIT COMPLEX (A - Z)
C*****precision > double
      DOUBLE PRECISION TEMP, TIMER
C*****END precision > double
C*****precision > single
C      REAL TEMP, TIMER
C*****END precision > single

      CALL CPUTIM (TEMP)
      TIMER = TEMP - TIMER

      RETURN
      END
      SUBROUTINE EXTENT (LENGTH, STRING)

C///////////////////////////////////////////////////////////////////////
C
C     EXTENT
C
C     FIND THE LAST NONBLANK CHARACTER IN A STRING.
C
C///////////////////////////////////////////////////////////////////////

      IMPLICIT COMPLEX (A - P, R - Z), INTEGER (Q)
      CHARACTER STRING*(*)
      INTEGER J, LENGTH

      LENGTH = 1
      DO 0100 J = 1, LEN (STRING)
         IF (STRING(J : J) .NE. ' ') LENGTH = J
0100  CONTINUE

      RETURN
      END
      SUBROUTINE GRAB (ERROR, TEXT, LAST, FIRST, NUMBER)

C///////////////////////////////////////////////////////////////////////
C
C     GRAB
C
C     RESERVE SPACE IN A DATA SPACE.
C
C///////////////////////////////////////////////////////////////////////

      IMPLICIT COMPLEX (A - Z)
      CHARACTER ID*9
      INTEGER FIRST, LAST, NUMBER, TEXT
      LOGICAL ERROR

      PARAMETER (ID = '  GRAB:  ')

C///  CHECK.

      ERROR = .NOT. (0 .LE. LAST)
      IF (ERROR) GO TO 9001

      ERROR = .NOT. (0 .LE. NUMBER)
      IF (ERROR) GO TO 9002

C///  GRAB THE SPACE.

      FIRST = LAST + 1
      LAST = LAST + MAX (1, NUMBER)

C///  ERROR MESSAGES.

      GO TO 99999

9001  IF (0 .LT. TEXT) WRITE (TEXT, 99001) ID, LAST
      GO TO 99999

9002  IF (0 .LT. TEXT) WRITE (TEXT, 99002) ID, NUMBER
      GO TO 99999

99001 FORMAT
     +   (/1X, A9, 'ERROR.  THE POINTERS FOR THE DATA SPACE ARE'
     +   /10X, 'INCONSISTENT.'
     +   //10X, I10, '  LAST ENTRY CURRENTLY IN USE')

99002 FORMAT
     +   (/1X, A9, 'ERROR.  THE NUMBER OF WORDS TO BE RESERVED IN THE'
     +   /10X, 'DATA SPACE MUST NOT BE NEGATIVE.'
     +   //10X, I10, '  NUMBER')

C///  EXIT.

99999 CONTINUE
      RETURN
      END
      SUBROUTINE LOGSTR (FORMAT, LENGTH, STRING, VALUE)

      IMPLICIT COMPLEX (A - Z)
      CHARACTER FORMAT*(*), STRING*(*), TEMP*80
      INTEGER FIRST, J, LAST, LENGTH
      LOGICAL FOUND
C*****precision > double
      DOUBLE PRECISION VALUE, ZERO
C*****END precision > double
C*****precision > single
C      REAL VALUE, ZERO
C*****END precision > single

C///  FORM THE WORD.

      ZERO = 0.0
      IF (ZERO .LT. VALUE) THEN
         WRITE (TEMP, FORMAT, ERR = 0200) LOG10 (VALUE)
      ELSE IF (VALUE .EQ. ZERO) THEN
         TEMP = 'NULL'
      ELSE
         TEMP = '?'
      END IF

C///  FIND THE FIRST AND LAST CHARACTERS.

      FOUND = .FALSE.
      LAST = 0
      DO 0100 J = 1, LEN (TEMP)
         IF (TEMP (J : J) .NE. ' ') THEN
            LAST = J
            IF (.NOT. FOUND) FIRST = J
            FOUND = .TRUE.
         END IF
0100  CONTINUE

C///  PACK THE WORD INTO THE STRING.

      IF (0 .LT. LENGTH .AND. LENGTH .LE. LEN (STRING)) THEN
         STRING (1 : LENGTH) = ' '
         IF (FOUND) THEN
            IF ((LAST - FIRST) .LT. LENGTH) THEN
               STRING (LENGTH - (LAST - FIRST) : LENGTH)
     +            = TEMP (FIRST : LAST)
            ELSE
               GO TO 0200
            END IF
         END IF
      ELSE
         GO TO 0200
      END IF

C///  SET THE STRING TO '!!! ...'.

      GO TO 99999

0200  CONTINUE
      DO 0300 J = 1, MIN (LEN (STRING), LENGTH)
         STRING (J : J) = '!'
0300  CONTINUE

C///  EXIT.

99999 CONTINUE
      RETURN
      END
      SUBROUTINE NEWTON
     +   (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
     +   EXIST, FUNCT, JACOB, LEVEL1, LEVEL2, N, REENT, REPORT, RETIRE,
     +   S, SHOW, SOLVE, STEPS, STMAX, STRIAL, SUCCES, X, XTRIAL, Y,
     +   YNORM, YTRIAL)

      IMPLICIT COMPLEX (A - P, R - Z), INTEGER (Q)
      CHARACTER COLUMN*16, HEADER*80, ID*9
      EXTERNAL COPY, LOGSTR, NORM
      INTEGER
     +   AGE, COUNT, ENTRY, FORMED, INDEX, J, K, LEVEL1, LEVEL2, LIMIT,
     +   N, REDUCE, REPORT, RETIRE, RETURN, ROUTE, STEPS, STMAX, TEXT
      LOGICAL
     +   ACCEPT, DECIDE, ERROR, EXIST, FORCE, FUNCT, JACOB, PRINTJ,
     +   PRINTS, PRINTY, REENT, SHOW, SOLVE, SUCCES
C*****precision > double
      DOUBLE PRECISION
C*****END precision > double
C*****precision > single
C      REAL
C*****END precision > single
     +   ABOVE, BELOW, BUFFER, COEFF, CONDIT, HALF, ONE, S, SNORM,
     +   STNORM, STRIAL, TEMP, VALUE, X, XTRIAL, Y, YNORM, YTNORM,
     +   YTRIAL, ZERO

      PARAMETER (ID = 'NEWTON:  ')
      PARAMETER (LIMIT = 5)

C     REPORT CODES
      PARAMETER
     +   (QNULL = 0, QDONE = 1, QDVRG = 2, QBNDS = 3, QSTDY = 4,
     +   QSMAX = 5)

      DIMENSION
     +   ABOVE(N), BELOW(N), BUFFER(N), COLUMN(7), HEADER(4, 2), S(N),
     +   STRIAL(N), X(N), XTRIAL(N), Y(N), YTRIAL(N)

C///  SAVE LOCAL VALUES DURING RETURNS FOR REVERSE COMMUNCIATION.

      SAVE

C     FORMED: SOLUTION INDEX AT WHICH A JACOBIAN WAS LAST COMPUTED.
C     INDEX: NUMBER OF THE LATEST ACCEPTED SOLUTION.

C///////////////////////////////////////////////////////////////////////
C
C     INITIALIZATION.
C
C///////////////////////////////////////////////////////////////////////

C///  EVERY-TIME INITIALIZATION.

C     TURN OFF ALL REVERSE COMMUNICATION FLAGS.
      DECIDE = .FALSE.
      FUNCT = .FALSE.
      JACOB = .FALSE.
      SHOW = .FALSE.
      SOLVE = .FALSE.

C     TURN OFF ALL COMPLETION STATUS FLAGS.
      ERROR = .FALSE.
      REPORT = QNULL
      SUCCES = .FALSE.

C     PRECISION-INDEPENDENT CONSTANTS.
      ONE = 1.0
      HALF = 0.5
      ZERO = 0.0

C///  IF THIS IS A RETURN CALL, THEN CONTINUE WHERE THE PROGRAM PAUSED.

      IF (REENT) THEN
         REENT = .FALSE.
         GO TO ROUTE
      END IF

C///  ONE-TIME INITIALIZATION.

      FORMED = - 1
      INDEX = 0
      PRINTJ = .FALSE.
      PRINTS = .FALSE.
      PRINTY = .FALSE.

C///  CHECK THE ARGUMENTS.

C     ERROR = .NOT. (1 .LE. N)
      IF (.NOT. (1 .LE. N)) ERROR = .TRUE.
      IF (ERROR) GO TO 9001

      COUNT = 0
      DO 1100 J = 1, N
         IF (.NOT. BELOW(J) .LT. ABOVE(J)) COUNT = COUNT + 1
1100  CONTINUE
C     ERROR = COUNT .NE. 0
      IF (COUNT .NE. 0) ERROR = .TRUE.
      IF (ERROR) GO TO 9002

      COUNT = 0
      DO 1200 J = 1, N
         IF (.NOT. (BELOW(J) .LE. X(J) .AND. X(J) .LE. ABOVE(J)))
     +      COUNT = COUNT + 1
1200  CONTINUE
C     ERROR = COUNT .NE. 0
      IF (COUNT .NE. 0) ERROR = .TRUE.
      IF (ERROR) GO TO 9003

C///  PRINT THE HEADER.

C                      123456789_123456789_123456789_123456789_12345
C                      123456   123456   123456   123456   123456
      HEADER (1, 1) = '         LOG10 MAX NORM                      '
      HEADER (2, 1) = '         ---------------------------------   '
      HEADER (3, 1) = 'SOLUTN   FUNCTN     STEP    TRIAL    TRIAL   '
      HEADER (4, 1) = 'NUMBER    VALUE     SIZE    VALUE     STEP   '

C                      123456789_12345
C                      123456   123456
      HEADER (1, 2) = '          LOG10'
      HEADER (2, 2) = ' LOG10   JACOBN'
      HEADER (3, 2) = 'DAMPNG   CONDTN'
      HEADER (4, 2) = ' COEFF   NUMBER'

      IF (LEVEL1 .GE. 1) THEN
         IF (TEXT .GT. 0) WRITE (TEXT, 10001)
     +      ID, ((HEADER(J, K), K = 1, 2), J = 1, 4)
      END IF

C///  EVALUATE THE FUNCTION AT THE INITIAL X AND STORE IN Y.

      PRINTY = .TRUE.
      ASSIGN 1300 TO RETURN
      GO TO 9911
1300  CONTINUE
      CALL NORM (N, YNORM, Y)

C///////////////////////////////////////////////////////////////////////
C
C     BLOCK FOR ONE STEP OF NEWTON'S ALGORITHM.
C
C///////////////////////////////////////////////////////////////////////

C///  FORM THE JACOBIAN AT THE LATEST X, AND RE-EVALUATE THE FUNCTION
C///  IN CASE THE FUNCTION HAS CHANGED, UNLESS THE ALGORITHM HAS JUST
C///  BEGUN AND A JACOBIAN ALREADY EXISTS.

      IF (INDEX .EQ. 0 .AND. EXIST) GO TO 2300

0200  CONTINUE
      FORMED = INDEX
      PRINTJ = .TRUE.
      ASSIGN 2100 TO RETURN
      GO TO 9921
2100  CONTINUE

      PRINTY = .TRUE.
      ASSIGN 2200 TO RETURN
      GO TO 9911
2200  CONTINUE
      CALL NORM (N, YNORM, Y)

2300  CONTINUE

C///  APPLY THE INVERSE JACOBIAN TO Y AND STORE THE STEP DIRECTION IN S.

      PRINTS = .TRUE.
      ASSIGN 2400 TO RETURN
      GO TO 9931
2400  CONTINUE
      CALL NORM (N, SNORM, S)

C///  IF THE SOLUTION IS ACCEPTABLE, THEN EXIT WITH SUCCESS.

      ASSIGN 2500 TO RETURN
      GO TO 9941
2500  CONTINUE
      ACCEPT = ACCEPT .AND. (1 .LE. INDEX)

      IF (.NOT. ACCEPT) GO TO 2800
         IF (LEVEL1 .GE. 1) THEN

C           PRINT A LINE FOR THE FINAL SOLUTION.
            DO 2600 J = 1, 7
               COLUMN(J) = ' '
2600        CONTINUE

            WRITE (COLUMN(1), '(I6)') INDEX
            CALL LOGSTR ('(F6.2)', 6, COLUMN(2), YNORM)
            CALL LOGSTR ('(F6.2)', 6, COLUMN(3), SNORM)
            IF (PRINTJ) CALL LOGSTR ('(F6.2)', 6, COLUMN(7), CONDIT)

            IF (TEXT .GT. 0) WRITE (TEXT, 10002) COLUMN

C           PRINT THE FINAL SOLUTION.
            IF (LEVEL1 .GE. LEVEL2 .AND. LEVEL2 .GE. 1) THEN
               IF (TEXT .GT. 0) WRITE (TEXT, 10003) ID
               ASSIGN 2700 TO RETURN
               IF (TEXT .GT. 0) GO TO 9951
            ELSE IF (LEVEL1 .GE. 1) THEN
               IF (TEXT .GT. 0) WRITE (TEXT, 10004) ID
            END IF
         END IF
2700     CONTINUE
         REPORT = QDONE
         SUCCES = .TRUE.
         GO TO 99999
2800  CONTINUE

C///  FIND THE LARGEST DAMPING COEFFICIENT BETWEEN 0 AND 1 THAT KEEPS
C///  THE TRIAL X WITHIN BOUNDS.  IF THE NEXT STEP WILL PLACE SOME
C///  VARIABLES AT THEIR UPPER OR LOWER BOUNDS, THEN MAKE PROVISIONS TO
C///  FORCE ONE OF THESE VARIABLES TO THE APPROPRIATE BOUNDARY.

0300  CONTINUE

      COEFF = ONE
      FORCE = .FALSE.
      DO 3100 J = 1, N
         IF (S(J) .GT. ZERO) THEN
            TEMP = (X(J) - BELOW(J)) / S(J)
            IF (TEMP .LE. COEFF) THEN
               COEFF = TEMP
               ENTRY = J
               FORCE = .TRUE.
               VALUE = BELOW(J)
            END IF
         ELSE IF (S(J) .LT. ZERO) THEN
            TEMP = (X(J) - ABOVE(J)) / S(J)
            IF (TEMP .LE. COEFF) THEN
               COEFF = TEMP
               ENTRY = J
               FORCE = .TRUE.
               VALUE = ABOVE(J)
            END IF
         END IF
3100  CONTINUE

C     ERROR = COEFF .LT. ZERO
      IF (COEFF .LT. ZERO) ERROR = .TRUE.
      IF (ERROR) GO TO 9004

      REDUCE = 0

C///  IF THE COEFFICIENT IS ZERO, THEN SOME VARIABLES LIE AT THEIR
C///  UPPER OR LOWER BOUNDS AND THE PRESENT STEP DIRECTION S WOULD
C///  CARRY THEM OUT OF BOUNDS.  S MUST BE BAD.  PRINT A LINE FOR THE
C///  FAILED S.

      IF (COEFF .EQ. ZERO) THEN

         IF (LEVEL1 .GE. 1) THEN
            DO 3200 J = 1, 7
               COLUMN(J) = ' '
3200        CONTINUE

            IF (PRINTY) WRITE (COLUMN(1), '(I6)') INDEX
            IF (PRINTY) CALL LOGSTR ('(F6.2)', 6, COLUMN(2), YNORM)
            IF (PRINTS) CALL LOGSTR ('(F6.2)', 6, COLUMN(3), SNORM)
            COLUMN(6) = '  ZERO'
            IF (PRINTJ) CALL LOGSTR ('(F6.2)', 6, COLUMN(7), CONDIT)

            IF (TEXT .GT. 0) WRITE (TEXT, 10002) COLUMN

            PRINTJ = .FALSE.
            PRINTS = .FALSE.
            PRINTY = .FALSE.
         END IF

C///  IF THE JACOBIAN IS NOT CURRENT, THEN CONTINUE THE ALGORITHM FROM
C///  THE POINT OF EVALUATING A NEW JACOBIAN AND FORMING A NEW S.
C///  OTHERWISE, THE ALGORITHM FAILS.

         IF (FORMED .LT. INDEX) GO TO 0200
         REPORT = QBNDS
         SUCCES = .FALSE.
         IF (LEVEL1 .GE. 1) THEN
            IF (TEXT .GT. 0) WRITE (TEXT, 10005) ID
            DO 0400 J = 1, N
               IF (BELOW(J) .EQ. X(J) .OR. X(J) .EQ. ABOVE(J)) THEN
                  IF (TEXT .GT. 0)
     +               WRITE (TEXT, 10006) J, BELOW(J), X(J), ABOVE(J)
               END IF
0400        CONTINUE
         END IF
         GO TO 99999
      END IF

C///  FORM THE TRIAL X, TAKING CARE THAT ROUNDING ERRORS DO NOT PLACE IT
C///  OUT OF BOUNDS.  BUT IF ONE OF THE VARIABLES SHOULD LIE AT A
C///  BOUNDARY, THEN THEN PLACE IT THERE BECAUSE ROUNDING ERRORS MAY
C///  LEAVE THE VARIABLE SHORT OF THE BOUNDARY.

0500  CONTINUE

      DO 5100 J = 1, N
         XTRIAL(J) = X(J) - COEFF * S(J)
5100  CONTINUE

      IF (REDUCE .EQ. 0) THEN
         DO 5200 J = 1, N
            XTRIAL(J) = MIN (XTRIAL(J), ABOVE(J))
5200     CONTINUE
         DO 5300 J = 1, N
            XTRIAL(J) = MAX (XTRIAL(J), BELOW(J))
5300     CONTINUE
      END IF

      IF (FORCE .AND. REDUCE .EQ. 0) THEN
         FORCE = .FALSE.
         XTRIAL(ENTRY) = VALUE
      END IF

C///  EVALUATE THE FUNCTION AT THE TRIAL X AND STORE IN THE TRIAL Y.

      ASSIGN 5400 TO RETURN
      GO TO 9961
5400  CONTINUE
      CALL NORM (N, YTNORM, YTRIAL)

C///  APPLY THE INVERSE JACOBIAN TO THE TRIAL Y AND STORE IN THE TRIAL
C///  S.

      ASSIGN 5500 TO RETURN
      GO TO 9971
5500  CONTINUE
      CALL NORM (N, STNORM, STRIAL)

C///////////////////////////////////////////////////////////////////////
C
C     PRINT.
C
C///////////////////////////////////////////////////////////////////////

      IF (LEVEL1 .GE. 1) THEN
         DO 6100 J = 1, 7
            COLUMN(J) = ' '
6100     CONTINUE

         IF (PRINTY) WRITE (COLUMN(1), '(I6)') INDEX
         IF (PRINTY) CALL LOGSTR ('(F6.2)', 6, COLUMN(2), YNORM)
         IF (PRINTS) CALL LOGSTR ('(F6.2)', 6, COLUMN(3), SNORM)
         CALL LOGSTR ('(F6.2)', 6, COLUMN(4), YTNORM)
         CALL LOGSTR ('(F6.2)', 6, COLUMN(5), STNORM)
         IF (COEFF .LT. ONE) CALL LOGSTR ('(F6.2)', 6, COLUMN(6), COEFF)
         IF (PRINTJ) CALL LOGSTR ('(F6.2)', 6, COLUMN(7), CONDIT)

         IF (TEXT .GT. 0) WRITE (TEXT, 10002) COLUMN

         PRINTJ = .FALSE.
         PRINTS = .FALSE.
         PRINTY = .FALSE.
      END IF

C///////////////////////////////////////////////////////////////////////
C
C     BLOCK FOR THE DECISION TREE.
C
C///////////////////////////////////////////////////////////////////////

C///  IF THE NORM OF THE TRIAL S IS SMALLER THAN THE NORM OF THE PRESENT
C///  S, THEN THE STEP IS SUCCESSFUL.  REPLACE X, Y, AND S BY THE TRIAL
C///  VALUES.

C     THIS IF-BLOCK IS IMPLEMENTED USING GO-TO'S IN ORDER TO PERMIT
C     THE CODE INSIDE THE BLOCK TO USE REVERSE COMMUNICATION.
C     IF (STNORM .LT. SNORM) THEN
      IF (.NOT. (STNORM .LT. SNORM)) GO TO 7500
         INDEX = INDEX + 1
         CALL COPY (N, XTRIAL, X)
         CALL COPY (N, YTRIAL, Y)
         CALL COPY (N, STRIAL, S)
         SNORM = STNORM
         YNORM = YTNORM
         PRINTS = .TRUE.
         PRINTY = .TRUE.

C///  IF THE SOLUTION IS ACCEPTABLE, THEN EXIT WITH SUCCESS.

         ASSIGN 7100 TO RETURN
         GO TO 9941
7100     CONTINUE
         ACCEPT = ACCEPT .AND. (1 .LE. INDEX)

         IF (.NOT. ACCEPT) GO TO 7400
            IF (LEVEL1 .GE. 1) THEN

C              PRINT A LINE FOR THE FINAL SOLUTION.
               DO 7200 J = 1, 7
                  COLUMN(J) = ' '
7200           CONTINUE

               WRITE (COLUMN(1), '(I6)') INDEX
               CALL LOGSTR ('(F6.2)', 6, COLUMN(2), YNORM)
               CALL LOGSTR ('(F6.2)', 6, COLUMN(3), SNORM)
               IF (PRINTJ) CALL LOGSTR ('(F6.2)', 6, COLUMN(7), CONDIT)

               IF (TEXT .GT. 0) WRITE (TEXT, 10002) COLUMN

C              PRINT THE FINAL SOLUTION.
               IF (LEVEL1 .GE. LEVEL2 .AND. LEVEL2 .GE. 1) THEN
                  IF (TEXT .GT. 0) WRITE (TEXT, 10003) ID
                  ASSIGN 7300 TO RETURN
                  IF (TEXT .GT. 0) GO TO 9951
               ELSE IF (LEVEL1 .GE. 1) THEN
                  IF (TEXT .GT. 0) WRITE (TEXT, 10004) ID
               END IF
            END IF
7300        CONTINUE
            REPORT = QDONE
            SUCCES = .TRUE.
            GO TO 99999
7400     CONTINUE

C///  OTHERWISE, IF THE LIMITING NUMBER OF STEPS HAVE BEEN TAKEN, THE
C///  ALGORITHM FAILS.

         IF (INDEX .EQ. STMAX) THEN
            REPORT = QSMAX
            SUCCES = .FALSE.
            IF (LEVEL1 .GE. 1) THEN
               IF (TEXT .GT. 0) WRITE (TEXT, 10007) ID
            END IF
            GO TO 99999
         END IF

C///  OTHERWISE, CONTINUE THE ALGORITHM.  IF THE JACOBIAN IS YOUNG,
C///  THEN CONTINUE AT THE POINT OF SELECTING THE NEW DAMPING
C///  COEFFICIENT.  BUT IF THE JACOBIAN IS OF RETIREMENT AGE, THEN
C///  CONTINUE AT THE POINT OF EVALUATING A NEW JACOBIAN.

         AGE = INDEX - MAX (0, FORMED)
         IF (AGE .LT. RETIRE) GO TO 0300
         GO TO 0200

C///  IF THE NORM OF THE TRIAL S IS NO SMALLER THAN THE NORM OF THE
C///  PRESENT S, THEN THE STEP IS NOT SUCCESSFUL.

C     ELSE
7500  CONTINUE

C///  IF THE DAMPING COEFFICIENT HAS BEEN REDUCED ONLY A FEW TIMES, THEN
C///  DECREASE THE COEFFICIENT AND CONTINUE THE ALGORITHM AT THE POINT
C///  OF FORMING A NEW TRIAL X.

         COEFF = HALF * COEFF
         REDUCE = REDUCE + 1
         IF (REDUCE .LE. LIMIT) GO TO 0500

C///  OTHERWISE, S MUST BE BAD.  IF THE JACOBIAN IS NOT CURRENT, THEN
C///  CONTINUE THE ALGORITHM FROM THE POINT OF EVALUATING A NEW
C///  JACOBIAN AND FORMING A NEW S.

         IF (FORMED .LT. INDEX) GO TO 0200

C///  OTHERWISE, THE ALGORITHM FAILS.

         REPORT = QDVRG
         SUCCES = .FALSE.
         IF (LEVEL1 .GE. 1) THEN
            IF (TEXT .GT. 0) WRITE (TEXT, 10008) ID
         END IF
         GO TO 99999
C     END IF

C///////////////////////////////////////////////////////////////////////
C
C     REQUEST REVERSE COMMUNICATION.
C
C///////////////////////////////////////////////////////////////////////

C///  EVALUATE THE FUNCTION AT X AND STORE IN Y.

9911  CONTINUE
      CALL COPY (N, X, BUFFER)
      FUNCT = .TRUE.
      REENT = .TRUE.
      ASSIGN 9912 TO ROUTE
      GO TO 99999
9912  CONTINUE
      CALL COPY (N, BUFFER, Y)
      GO TO RETURN

C///  FORM THE JACOBIAN AT THE LATEST X.

9921  CONTINUE
      CALL COPY (N, X, BUFFER)
      JACOB = .TRUE.
      REENT = .TRUE.
      ASSIGN 9922 TO ROUTE
      GO TO 99999
9922  CONTINUE
      GO TO RETURN

C///  APPLY THE INVERSE JACOBIAN TO Y AND STORE THE STEP DIRECTION IN S.

9931  CONTINUE
      CALL COPY (N, Y, BUFFER)
      SOLVE = .TRUE.
      REENT = .TRUE.
      ASSIGN 9932 TO ROUTE
      GO TO 99999
9932  CONTINUE
      CALL COPY (N, BUFFER, S)
      GO TO RETURN

C///  CHECK WHETHER THE NEW S IS SUFFICIENTLY SMALL THAT THE SOLUTION IS
C///  ACCEPTABLE.

9941  CONTINUE
      DECIDE = .TRUE.
      REENT = .TRUE.
      ASSIGN 9942 TO ROUTE
      GO TO 99999
9942  CONTINUE
      GO TO RETURN

C///  PRINT THE SOLUTION.

9951  CONTINUE
      CALL COPY (N, X, BUFFER)
      SHOW = .TRUE.
      REENT = .TRUE.
      ASSIGN 9952 TO ROUTE
      GO TO 99999
9952  CONTINUE
      GO TO RETURN

C///  EVALUATE THE FUNCTION AT THE TRIAL X AND STORE IN THE TRIAL Y.

9961  CONTINUE
      CALL COPY (N, XTRIAL, BUFFER)
      FUNCT = .TRUE.
      REENT = .TRUE.
      ASSIGN 9962 TO ROUTE
      GO TO 99999
9962  CONTINUE
      CALL COPY (N, BUFFER, YTRIAL)
      GO TO RETURN

C///  APPLY THE INVERSE JACOBIAN TO THE TRIAL Y AND STORE THE TRIAL
C///  STEP DIRECTION IN THE TRIAL S.

9971  CONTINUE
      CALL COPY (N, YTRIAL, BUFFER)
      SOLVE = .TRUE.
      REENT = .TRUE.
      ASSIGN 9972 TO ROUTE
      GO TO 99999
9972  CONTINUE
      CALL COPY (N, BUFFER, STRIAL)
      GO TO RETURN

C///////////////////////////////////////////////////////////////////////
C
C     FORMAT STATEMENTS FOR INFORMATIVE MESSAGES.
C
C///////////////////////////////////////////////////////////////////////

10001 FORMAT
     +   (/1X, A9, 'LOG OF THE DAMPED AND SIMPLIFIED NEWTON ALGORITHM.'
     +   /4(/10X, A45, A15)/)

10002 FORMAT
     +   (10X, A6, 6(3X, A6))

10003 FORMAT
     +   (/1X, A9, 'SUCCESS.  THE SOLUTION:')

10004 FORMAT
     +   (/1X, A9, 'SUCCESS.')

10005 FORMAT
     +   (/1X, A9, 'FAILURE.  SOME VARIABLES LIE AT THEIR UPPER OR'
     +   /10X, 'LOWER BOUNDS.  THE NEXT STEP WOULD CARRY THEM OUT OF'
     +   /10X, 'BOUNDS.'
C               123456789   123456789   123456789   123456789
     +  //10X, ' VARIABLE       LOWER                   UPPER'
     +   /10X, '   NUMBER       BOUND       VALUE       BOUND'
     +   /)

10006 FORMAT
     +   (10X, I9, 3(3X, E9.2))

10007 FORMAT
     +   (/1X, A9, 'FAILURE.  THE LIMITING NUMBER OF STEPS HAS BEEN'
     +   /10X, 'TAKEN.')

10008 FORMAT
     +   (/1X, A9, 'FAILURE.  NO DAMPING COEFFICIENT PRODUCED A TRIAL'
     +   /10X, 'SOLUTION WITH A STEP VECTOR SHORTER THAN THE PRESENT'
     +   /10X, 'UNDAMPED STEP.')

C///////////////////////////////////////////////////////////////////////
C
C     ERROR MESSAGES.
C
C///////////////////////////////////////////////////////////////////////

9001  IF (TEXT .GT. 0)
     +   WRITE (TEXT, 99001) ID, N
      GO TO 99999

9002  IF (TEXT .GT. 0) THEN
         WRITE (TEXT, 99002) ID, N, COUNT
         COUNT = 0
         DO 9100 J = 1, N
            IF (.NOT. BELOW(J) .LT. ABOVE(J)) THEN
               COUNT = COUNT + 1
               IF (COUNT .LE. 20) THEN
                  WRITE (TEXT, 98001) J, BELOW(J), ABOVE(J)
               ELSE IF (COUNT .EQ. 21) THEN
                  WRITE (TEXT, 98002)
               END IF
            END IF
9100     CONTINUE
      END IF
      GO TO 99999

9003  IF (TEXT .GT. 0)
     +   WRITE (TEXT, 99003) ID, COUNT
      GO TO 99999

9004  IF (TEXT .GT. 0)
     +   WRITE (TEXT, 99004) ID, INDEX, COEFF
      GO TO 99999

99001 FORMAT
     + (/1X, A9, 'ERROR.  THE NUMBER OF VARIABLES IS NOT POSITIVE.'
     + //1X, I10, '  NUMBER')

99002 FORMAT
     + (/1X, A9, 'ERROR.  THE LEFT BOUNDS UPON SOME SOLUTION VALUES'
     +  /10X, 'EQUAL OR EXCEED THE RIGHT BOUNDS.'
     + //1X, I10, '  SOLUTION VARIABLES'
     +  /1X, I10, '  VARIABLES WITH INCONSISENT BOUNDS'
C             123456789_  123456789_  123456789_
     + //1X, '     INDEX  LEFT BOUND       RIGHT')

98001 FORMAT
     + (1X, I10, 1P, 2X, E10.2, 2X, E10.2)

98002 FORMAT
     + (1X, '       ...')

99003 FORMAT
     + (/1X, A9, 'ERROR.  SOME INITIAL SOLUTION VALUES LIE OUTSIDE'
     +  /10X, 'THEIR BOUNDS.'
     + //1X, I10, '  VALUES OUT OF BOUNDS')

99004 FORMAT
     + (/1X, A9, 'ERROR.  THE INITIAL DAMPING COEFFICIENT IS NEGATIVE.'
     +  /10X, 'THE PROGRAM LOGIC SHOULD NOT PERMIT THIS.  EITHER THE'
     +  /10X, 'LOGIC IS WRONG, OR THE MACHINE ARITHMETIC IS'
     +  /10X, 'INCONSISTENT.'
     + //1X, I10, '  LATEST SOLUTION VECTOR'
     +  /1X, E10.2, '  COEFFICIENT')

C///  EXIT.

99999 CONTINUE

C     COPY THE PROTECTED LOCAL VARIABLE INTO THE ARGUMENT LIST.
      STEPS = INDEX

      RETURN
      END
      SUBROUTINE NORM (N, VALUE, X)

      IMPLICIT COMPLEX (A - Z)
C*****precision > double
      DOUBLE PRECISION VALUE, X
C*****END precision > double
      INTEGER J, N
C*****precision > single
C      REAL VALUE, X
C*****END precision > single

      DIMENSION X(N)

      VALUE = 0.0
      DO 0100 J = 1, N
         VALUE = MAX (VALUE, ABS (X(J)))
0100  CONTINUE

      RETURN
      END
      SUBROUTINE REFINE
     +   (ERROR, TEXT, ABOVE, ACTIVE, BELOW, BUFFER, COMPS, LEVEL1,
     +   LEVEL2, MARK, MESH, NEWMSH, PLIMIT, PMAX, POINTS, RATIO1,
     +   RATIO2, REENT, SHOW, SUCCES, TOLER0, TOLER1, TOLER2, UPDATE,
     +   VARY, VARY1, VARY2, X)

      IMPLICIT COMPLEX (A - Z)
      CHARACTER ID*9
      EXTERNAL COPY
      INTEGER
     +   ACT, ACTSIG, BOUND, COMPS, COUNT, COUNT1, COUNT2, EQUAL,
     +   FORMER, HIGHER, IDEAL, J, K, LEVEL1, LEVEL2, MORE, NEW, OLD,
     +   PLIMIT, PMAX, PMIN, POINTS, ROUTE, STRICT, TEXT, TOTAL, VARY,
     +   VARY1, VARY2
      LOGICAL
     +   ACTIVE, ERROR, MARK, NEWMSH, REENT, SHOW, SIGNIF, SUCCES,
     +   UPDATE
C*****precision > double
      DOUBLE PRECISION
     +   ABOVE, BELOW, BUFFER, DIFFER, HALF, LEFT, LENGTH, LOWER,
     +   MAXMAG, MEAN, MESH, ONE, RANGE, RATIO1, RATIO2, RIGHT, TEMP,
     +   TEMP1, TEMP2, TOLER0, TOLER1, TOLER2, UPPER, X, ZERO
C*****END precision > double
C*****precision > single
C      REAL
C     +   ABOVE, BELOW, BUFFER, DIFFER, HALF, LEFT, LENGTH, LOWER,
C     +   MAXMAG, MEAN, MESH, ONE, RANGE, RATIO1, RATIO2, RIGHT, TEMP,
C     +   TEMP1, TEMP2, TOLER0, TOLER1, TOLER2, UPPER, X, ZERO
C*****END precision > single

      PARAMETER (ID = 'REFINE:  ')
      PARAMETER (PMIN = 3)

      DIMENSION
     +   ABOVE(COMPS, PMAX), ACTIVE(COMPS), BELOW(COMPS, PMAX),
     +   BUFFER(COMPS * PMAX), MARK(PMAX),
     +   MESH(PMAX), RATIO1(PMAX), RATIO2(PMAX), VARY(PMAX),
     +   VARY1(PMAX), VARY2(PMAX), X(COMPS, PMAX)

C///  SAVE LOCAL VARIABLES DURING RETURNS FOR REVERSE COMMUNICATION.

      SAVE

C///////////////////////////////////////////////////////////////////////
C
C     RETURN FROM REVERSE COMMUNICATION.
C
C///////////////////////////////////////////////////////////////////////

C     TURN OFF ALL REVERSE COMMUNICATION FLAGS.
      SHOW = .FALSE.
      UPDATE = .FALSE.

C     TURN OFF ALL COMPLETION STATUS FLAGS.
      ERROR = .FALSE.
      NEWMSH = .FALSE.
      SUCCES = .FALSE.

C///  CONTINUE WHERE THE PROGRAM PAUSED.

      IF (REENT) THEN
         REENT = .FALSE.
         GO TO ROUTE
      END IF

C///////////////////////////////////////////////////////////////////////
C
C     ENTER.
C
C///////////////////////////////////////////////////////////////////////

C///  ONE-TIME INITIALIZATION.

C     PRECISION-INDEPENDENT CONSTANTS.
      HALF = 0.5
      ONE = 1.0
      ZERO = 0.0

C///  LEVEL1 PRINTING.

      IF (LEVEL1 .GE. 1 .AND. TEXT .GT. 0) WRITE (TEXT, 10001) ID

C///  CHECK THE ARGUMENTS.

      ERROR = .NOT. (0 .LE. PLIMIT)
      IF (ERROR) GO TO 9001

      ERROR = .NOT. (PMIN .LE. POINTS .AND. POINTS .LE. PMAX)
      IF (ERROR) GO TO 9002

      COUNT = 0
      DO 1100 J = 1, COMPS
         IF (ACTIVE(J)) COUNT = COUNT + 1
1100  CONTINUE
      ERROR = .NOT. (1 .LE. COUNT)
      IF (ERROR) GO TO 9003

      ERROR = .NOT. (ZERO .LE. TOLER0)
      IF (ERROR) GO TO 9004

      ERROR = .NOT. (ZERO .LE. TOLER1 .AND. TOLER1 .LE. ONE
     +         .AND. ZERO .LE. TOLER2 .AND. TOLER2 .LE. ONE)
      IF (ERROR) GO TO 9005

      COUNT = 0
      DO 1200 K = 1, POINTS - 1
         IF (MESH(K) .LT. MESH(K + 1)) COUNT = COUNT + 1
1200  CONTINUE
      ERROR = .NOT. (COUNT .EQ. 0 .OR. COUNT .EQ. POINTS - 1)
      IF (ERROR) GO TO 9006

      LENGTH = ABS (MESH(POINTS) - MESH(1))
      COUNT = 0
      DO 1300 K = 1, POINTS - 1
         DIFFER = ABS (MESH(K) - MESH(K + 1))
         TEMP = LENGTH + DIFFER
         IF (TEMP .EQ. LENGTH) COUNT = COUNT + 1
1300  CONTINUE
      ERROR = .NOT. (COUNT .EQ. 0)
      IF (ERROR) GO TO 9007

C///////////////////////////////////////////////////////////////////////
C
C     AT EACH INTERVAL, COUNT THE ACTIVE, SIGNIFICANT COMPONENTS THAT
C     VARY TOO GREATLY.
C
C///////////////////////////////////////////////////////////////////////

      ACT = 0
      ACTSIG = 0

      DO 2100 K = 1, POINTS
         MARK(K) = .FALSE.
         RATIO1(K) = ZERO
         RATIO2(K) = ZERO
         VARY1(K) = 0
         VARY2(K) = 0
2100  CONTINUE

C///  TOP OF THE LOOP OVER THE COMPONENTS.

      DO 2600 J = 1, COMPS
         IF (ACTIVE(J)) THEN
            ACT = ACT + 1

C///  FIND THE RANGE AND MAXIMUM MAGNITUDE OF THE COMPONENT.

            LOWER = X(J, 1)
            UPPER = X(J, 1)
            MAXMAG = ABS (X(J, 1))
            DO 2200 K = 2, POINTS
               LOWER = MIN (LOWER, X(J, K))
               UPPER = MAX (UPPER, X(J, K))
               MAXMAG = MAX (MAXMAG, ABS (X(J, K)))
2200        CONTINUE
            RANGE = UPPER - LOWER

C///  DECIDE WHETHER THE COMPONENT IS SIGNIFICANT.

            SIGNIF = ABS (RANGE) .GT. TOLER0 * MAX (ONE, MAXMAG)
            IF (SIGNIF) THEN
               ACTSIG = ACTSIG + 1

C///  AT EACH INTERVAL, SEE WHETHER THE COMPONENT'S CHANGE EXCEEDS SOME
C///  FRACTION OF THE COMPONENT'S GLOBAL CHANGE.

               DO 2300 K = 1, POINTS - 1
                  DIFFER = ABS (X(J, K + 1) - X(J, K))
                  IF (RANGE .NE. ZERO)
     +               RATIO1(K) = MAX (RATIO1(K), DIFFER / RANGE)
                  IF (DIFFER .GT. TOLER1 * RANGE)
     +               VARY1(K) = VARY1(K) + 1
2300           CONTINUE

C///  FIND THE GLOBAL CHANGE OF THE COMPONENT'S DERIVATIVE.

               TEMP = (X(J, 2) - X(J, 1)) / (MESH(2) - MESH(1))
               LOWER = TEMP
               UPPER = TEMP
               DO 2400 K = 2, POINTS - 1
                  TEMP = (X(J, K + 1) - X(J, K))
     +                 / (MESH(K + 1) - MESH(K))
                  LOWER = MIN (LOWER, TEMP)
                  UPPER = MAX (UPPER, TEMP)
2400           CONTINUE
               RANGE = UPPER - LOWER

C///  AT EACH INTERIOR POINT, SEE WHETHER THE DERIVATIVE'S CHANGE
C///  EXCEEDS SOME FRACTION OF THE DERIVATIVE'S GLOBAL CHANGE.

               DO 2500 K = 2, POINTS - 1
                  LEFT = (X(J, K) - X(J, K - 1))
     +                 / (MESH(K) - MESH(K - 1))
                  RIGHT = (X(J, K + 1) - X(J, K))
     +                  / (MESH(K + 1) - MESH(K))
                  DIFFER = ABS (LEFT - RIGHT)
                  IF (RANGE .NE. ZERO)
     +               RATIO2(K) = MAX (RATIO2(K), DIFFER / RANGE)
                  IF (DIFFER .GT. TOLER2 * RANGE)
     +               VARY2(K) = VARY2(K) + 1
2500           CONTINUE

C///  BOTTOM OF THE LOOP OVER THE COMPONENTS.

            END IF
         END IF
2600  CONTINUE

C///  COUNT THE INTERVALS IN WHICH VARIATIONS THAT ARE TOO LARGE OCCUR.

      IDEAL = 0
      DO 2700 K = 1, POINTS - 1
         VARY(K) = VARY1(K)
         IF (1 .LT. K) VARY(K) = VARY(K) + VARY2(K)
         IF (K .LT. POINTS - 1) VARY(K) = VARY(K) + VARY2(K + 1)
         IF (0 .LT. VARY(K)) IDEAL = IDEAL + 1
2700  CONTINUE

C///////////////////////////////////////////////////////////////////////
C
C     SELECT THE INTERVALS TO HALVE.
C
C     HALVE INTERVALS ON WHICH COMPONENTS OR DERIVATIVES VARY TOO
C     GREATLY.  IF IT IS NOT POSSIBLE TO HALVE ALL SUCH INTERVALS, THEN
C     FIND A LOWER BOUND FOR THE NUMBER OF VARYING COMPONENTS PER
C     INTERVAL SO THAT ALL INTERVALS WITH STRICTLY MORE CAN BE HALVED.
C
C///////////////////////////////////////////////////////////////////////

C///  FIND THE BOUND.

      BOUND = 0
      MORE = MAX (0, MIN (PMAX - POINTS, PLIMIT))

3100  CONTINUE
      EQUAL = 0
      STRICT = 0
      DO 3200 K = 1, POINTS - 1
         IF (VARY(K) .EQ. BOUND) THEN
            EQUAL = EQUAL + 1
         ELSE IF (VARY(K) .GT. BOUND) THEN
            STRICT = STRICT + 1
            IF (STRICT .EQ. 1) THEN
               HIGHER = VARY(K)
            ELSE
               HIGHER = MIN (HIGHER, VARY(K))
            END IF
         END IF
3200  CONTINUE

      IF (STRICT .GT. MORE) THEN
         BOUND = HIGHER
         GO TO 3100
      END IF

C///  DETERMINE HOW MANY INTERVALS TO HALVE OF THOSE WHOSE NUMBER OF
C///  VARIATIONS EXACTLY EQUAL THE BOUND.

C     IF (0 .LT. BOUND .AND. 0 .EQ. STRICT) THEN
C        EQUAL = MIN (EQUAL, MORE)
C     ELSE
C        EQUAL = 0
C     END IF

      IF (0 .LT. BOUND) THEN
         EQUAL = MIN (EQUAL, MORE - STRICT)
      ELSE
         EQUAL = 0
      END IF

C///  MARK THE INTERVALS TO HALVE.

      COUNT = 0
      DO 3300 K = 1, POINTS - 1
         IF (BOUND .LT. VARY(K)) THEN
            MARK(K) = .TRUE.
         ELSE IF (BOUND .EQ. VARY(K) .AND. COUNT .LT. EQUAL) THEN
            COUNT = COUNT + 1
            MARK(K) = .TRUE.
         ELSE
            MARK(K) = .FALSE.
         END IF
3300  CONTINUE

C///////////////////////////////////////////////////////////////////////
C
C     HALVE THE INTERVALS, IF ANY.
C
C///////////////////////////////////////////////////////////////////////

C///  FORM THE TOTAL NUMBER OF POINTS IN THE NEW MESH.

      TOTAL = POINTS + EQUAL + STRICT
      IF (0 .EQ. EQUAL + STRICT) GO TO 4700

C///  TOP OF THE BLOCK TO CREATE THE NEW MESH.  CHECK THE INTERVALS.

      COUNT1 = 0
      COUNT2 = 0
      LENGTH = ABS (MESH(POINTS) - MESH(1))
      DO 4100 K = 1, POINTS - 1
         IF (MARK(K)) THEN
            MEAN = HALF * (MESH(K) + MESH(K + 1))

            TEMP = LENGTH + ABS (MESH(K + 1) - MEAN)
            IF (TEMP .EQ. LENGTH) COUNT1 = COUNT1 + 1
            TEMP = LENGTH + ABS (MEAN - MESH(K))
            IF (TEMP .EQ. LENGTH) COUNT1 = COUNT1 + 1

            IF (.NOT. ((MESH(K) .LT. MEAN .AND. MEAN .LT. MESH(K + 1))
     +         .OR. (MESH(K + 1) .LT. MEAN .AND. MEAN .LT. MESH(K))))
     +         COUNT2 = COUNT2 + 1
         END IF
4100  CONTINUE
      ERROR = .NOT. (COUNT1 .EQ. 0 .AND. COUNT2 .EQ. 0)
      IF (ERROR) GO TO 9008

C///  ADD THE NEW POINTS.  INTERPOLATE X AND THE BOUNDS.

      NEW = TOTAL
      DO 4400 OLD = POINTS, 2, - 1
         MESH(NEW) = MESH(OLD)
         DO 4200 J = 1, COMPS
            ABOVE(J, NEW) = ABOVE(J, OLD)
            BELOW(J, NEW) = BELOW(J, OLD)
            X(J, NEW) = X(J, OLD)
4200     CONTINUE
         NEW = NEW - 1

         IF (MARK(OLD - 1)) THEN
            MESH(NEW) = HALF * (MESH(OLD) + MESH(OLD - 1))
            DO 4300 J = 1, COMPS
               ABOVE(J, NEW)
     +            = HALF * (ABOVE(J, OLD) + ABOVE(J, OLD - 1))
               BELOW(J, NEW)
     +            = HALF * (BELOW(J, OLD) + BELOW(J, OLD - 1))
               X(J, NEW) = HALF * (X(J, OLD) + X(J, OLD - 1))
4300        CONTINUE
            NEW = NEW - 1
         END IF
4400  CONTINUE

C///  MARK THE NEW POINTS.

      NEW = TOTAL
      DO 4500 OLD = POINTS, 2, - 1
         MARK(NEW) = .FALSE.
         NEW = NEW - 1
         IF (MARK(OLD - 1)) THEN
            MARK(NEW) = .TRUE.
            NEW = NEW - 1
         END IF
4500  CONTINUE
      MARK(NEW) = .FALSE.

C///  UPDATE THE NUMBER OF POINTS.

      FORMER = POINTS
      POINTS = TOTAL

C///  ALLOW THE USER TO UPDATE THE SOLUTION AND BOUNDS.

      CALL COPY (COMPS * POINTS, X, BUFFER)
      REENT = .TRUE.
      UPDATE = .TRUE.
      ASSIGN 4600 TO ROUTE
      GO TO 99999
4600  CONTINUE
      CALL COPY (COMPS * POINTS, BUFFER, X)

C///  BOTTOM OF THE BLOCK TO CREATE A NEW MESH.

4700  CONTINUE

C///////////////////////////////////////////////////////////////////////
C
C     PRINT.
C
C///////////////////////////////////////////////////////////////////////

C///  LEVEL1 PRINTING.  SHOW THE COUNTS AND TOLERANCES.

      IF (LEVEL1 .GE. 1 .AND. TEXT .GT. 0) THEN
         WRITE (TEXT, 10002) COMPS - ACT, ACT - ACTSIG

         IF (ACTSIG .EQ. 0) THEN
            WRITE (TEXT, 10003) ID
         ELSE
            TEMP1 = RATIO1(1)
            DO 5100 K = 2, FORMER - 1
               TEMP1 = MAX (TEMP1, RATIO1(K))
5100        CONTINUE

            TEMP2 = RATIO2(2)
            DO 5200 K = 3, FORMER - 1
               TEMP2 = MAX (TEMP2, RATIO2(K))
5200        CONTINUE

            WRITE (TEXT, 10004)
     +         TEMP1, TOLER1, TEMP2, TOLER2, IDEAL, EQUAL + STRICT

            IF (IDEAL .GT. 0 .AND. MORE .EQ. 0) WRITE (TEXT, 10005)

C///  PRINT THE MESH.

            IF (IDEAL .EQ. 0) THEN
               WRITE (TEXT, 10006)
            ELSE
               WRITE (TEXT, 10007)
            END IF

            WRITE (TEXT, 10008)

            OLD = 0
            DO 5300 K = 1, POINTS
               IF (.NOT. MARK(K)) THEN
                  OLD = OLD + 1
                  IF (1 .LT. K) THEN
                     IF (MARK(K - 1)) THEN
                        WRITE (TEXT, 10009)
     +                     K - 1, MESH(K - 1),
     +                     RATIO1(OLD - 1), VARY1(OLD - 1)
                     ELSE
                        WRITE (TEXT, 10010)
     +                     RATIO1(OLD - 1), VARY1(OLD - 1)
                     END IF
                  END IF
                  IF (1 .LT. K .AND. K .LT. POINTS) THEN
                     WRITE (TEXT, 10011)
     +                  K, MESH(K), RATIO2(OLD), VARY2(OLD)
                  ELSE
                     WRITE (TEXT, 10011) K, MESH(K)
                  END IF
               END IF
5300        CONTINUE

         END IF
      END IF

C///  LEVEL2 PRINTING.

      IF (1 .LE. LEVEL2 .AND. LEVEL2 .LE. LEVEL1 .AND.
     +   0 .LT. EQUAL + STRICT .AND. TEXT .GT. 0) THEN
         WRITE (TEXT, 10012) ID
         CALL COPY (COMPS * POINTS, X, BUFFER)
         SHOW = .TRUE.
         REENT = .TRUE.
         ASSIGN 5400 TO ROUTE
         GO TO 99999
      END IF
5400  CONTINUE

C///////////////////////////////////////////////////////////////////////
C
C     DEFINE COMPLETION STATUS FLAGS.
C
C///////////////////////////////////////////////////////////////////////

      SUCCES = 0 .EQ. IDEAL
      NEWMSH = 0 .LT. EQUAL + STRICT

C///////////////////////////////////////////////////////////////////////
C
C     INFORMATIVE MESSAGES.
C
C///////////////////////////////////////////////////////////////////////

10001 FORMAT
     +  (/1X, A9, 'REFINE THE MESH, IF NECESSARY.')

10002 FORMAT
     +  (/10X, I10, '  INACTIVE COMPONENTS'
     +   /10X, I10, '  INSIGNIFICANT, ACTIVE COMPONENTS')

10003 FORMAT
     +  (/1X, A9, 'WARNING.  THERE ARE NO ACTIVE, SIGNIFICANT'
     +   /10X, 'COMPONENTS.')

10004 FORMAT
     +  (/10X, F10.5, '  LARGEST CHANGE IN VALUE (LOCAL / GLOBAL)'
     +   /10X, F10.5, '  DESIRED BOUND'
     +  //10X, F10.5, '  LARGEST CHANGE IN DERIVATIVE (LOCAL / GLOBAL)'
     +   /10X, F10.5, '  DESIRED BOUND'
     +  //10X, I10, '  POTENTIAL NUMBER OF NEW MESH POINTS'
     +   /10X, I10, '  ACTUAL NUMBER')

10005 FORMAT
     +  (/10X, 'WARNING.  MORE POINTS WOULD BE ADDED TO THE MESH WERE'
     +   /10X, 'IT NOT FOR THE LIMIT ON NEW POINTS OR THE LIMIT ON'
     +   /10X, 'TOTAL POINTS.')

10006 FORMAT
     +  (/10X, 'NO NEW POINTS ARE NEEDED.  THE MESH:')

10007 FORMAT
     +  (/10X, 'THE MESH, "*" DENOTES A NEW POINT:')

10008 FORMAT
C                     123                123    1    123    1
C               123456   1234567890123456   1234 1234   1234 1234
     +  (/10X, '                              LARGEST   NUMBER OF'
     +   /10X, '                            VARIATION   TOO LARGE'
     +   /10X, ' INDEX         MESH POINT       RATIO      RATIOS'
     +   /10X, '------   ----------------   ---------   ---------'
     +   /10X, '                             1ST  2ND    1ST  2ND')

10009 FORMAT
     +  (10X, I6, '*  ', E16.9, 3X, F4.2, 8X, I4)

10010 FORMAT
     +  (38X, F4.2, 8X, I4)

10011 FORMAT
     +  (10X, I6, '   ', E16.9, 8X, F4.2, 8X, I4)

10012 FORMAT
     +  (/1X, A9, 'THE SOLUTION ESTIMATE ON THE NEW MESH:')

C///////////////////////////////////////////////////////////////////////
C
C     ERROR MESSAGES.
C
C///////////////////////////////////////////////////////////////////////

      GO TO 99999

9001  IF (TEXT .GT. 0) WRITE (TEXT, 99001)
     +   ID, PLIMIT
      GO TO 99999

9002  IF (TEXT .GT. 0) WRITE (TEXT, 99002)
     +   ID, POINTS, PMIN, PMAX
      GO TO 99999

9003  IF (TEXT .GT. 0) WRITE (TEXT, 99003)
     +   ID, COUNT, COMPS
      GO TO 99999

9004  IF (TEXT .GT. 0) WRITE (TEXT, 99004)
     +   ID, TOLER0
      GO TO 99999

9005  IF (TEXT .GT. 0) WRITE (TEXT, 99005)
     +   ID, TOLER1, TOLER2
      GO TO 99999

9006  IF (TEXT .GT. 0) WRITE (TEXT, 99006)
     +   ID, COUNT, POINTS - 1 - COUNT
      GO TO 99999

9007  IF (TEXT .GT. 0) WRITE (TEXT, 99007)
     +   ID, LENGTH, COUNT
      GO TO 99999

9008  IF (TEXT .GT. 0) WRITE (TEXT, 99008) ID, COUNT1, COUNT2
      GO TO 99999

99001 FORMAT
     +  (/1X, A9, 'ERROR.  THE LIMIT ON THE NUMBER OF NEW POINTS MUST'
     +   /10X, 'BE POSITIVE.'
     +   //10X, I10, '  LIMIT')

99002 FORMAT
     +  (/1X, A9, 'ERROR.  THERE ARE TOO MANY OR TOO FEW POINTS.'
     +  //10X, I10, '  POINTS'
     +   /10X, I10, '  MINIMUM (SET BY THIS SUBROUTINE)'
     +   /10X, I10, '  MAXIMUM')

99003 FORMAT
     +  (/1X, A9, 'ERROR.  THERE ARE NO ACTIVE COMPONENTS.'
     +  //10X, I10, '  ACTIVE COMPONENTS'
     +   /10X, I10, '  TOTAL COMPONENTS')

99004 FORMAT
     +  (/1X, A9, 'ERROR.  THE SIGNIFICANCE LEVEL IS NEGATIVE.'
     +  //10X, E10.2, '  TOLER0')

99005 FORMAT
     +  (/1X, A9, 'ERROR.  THE BOUNDS UPON THE RELATIVE VARIATIONS'
     +   /10X, 'DO NOT LIE BETWEEN ZERO AND ONE, INCLUSIVE.'
     +  //10X, E10.2, '  TOLER1'
     +   /10X, E10.2, '  TOLER2')

99006 FORMAT
     +  (/1X, A9, 'ERROR.  THE POINTS ARE NOT ORDERED.'
     +  //10X, I10, '  INTERVALS WITH ENDPOINTS ORDERED LEFT TO RIGHT'
     +   /10X, I10, '  INTERVALS WITH ENDPOINTS ORDERED RIGHT TO LEFT')

99007 FORMAT
     +  (/1X, A9, 'ERROR.  SOME INTERVALS ARE SMALLER THAN MACHINE'
     +   /10X, 'EPSILON TIMES THE TOTAL LENGTH OF ALL INTERVALS.'
     +  //10X, E10.2, '  TOTAL LENGTH'
     +   /10X, I10, '  SMALL INTERVALS')

99008 FORMAT
     +  (/1X, A9, 'ERROR.  EITHER (1) SOME NEW INTERVALS WOULD BE'
     +   /10X, 'SMALLER THAN MACHINE EPSILON TIMES THE TOTAL LENGTH OF'
     +   /10X, 'ALL INTERVALS, OR (2) ROUNDING ERRORS PLACE SOME NEW'
     +   /10X, 'POINTS OUTSIDE THE INTERVALS THEY SHOULD HALVE.'
     +  //10X, I10, '  NEW INTERVALS TOO SMALL'
     +   /10X, I10, '  NEW POINTS OUT OF ORDER')

99999 CONTINUE
      RETURN
      END
      SUBROUTINE TIMSTP
     +   (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
     +   DESIRE, FUNCT, JACOB, LEVEL1, LEVEL2, N, REENT, REPORT,
     +   RETIRE, S, SHOW, SIZE1, SOLVE, STEP, STMAX, STORE, STRIAL,
     +   SUCCES, TDAGE, TDEC, TIME, TINC, TMAX, TMIN, TOLER, X, XSAVE,
     +   XTRIAL, Y, YNORM, YTRIAL)

      IMPLICIT COMPLEX (A - P, R - Z), INTEGER (Q)
      CHARACTER COLUMN*16, HEADER*80, ID*9, REMARK*80, STRING*80
      EXTERNAL COPY, EXTENT, LOGSTR, NEWTON, NORM
      INTEGER
     +   AGE, BEGIN, COUNT, DESIRE, J, K, LENGTH, LEVEL1, LEVEL2, N,
     +   NUMBER, REPORT, RETIRE, ROUTE, STEP, STMAX, TDAGE, TEXT,
     +   TRIALS, TRMAX, XREPOR
      LOGICAL
     +   ACCEPT, AGAIN, DECIDE, ERROR, EXIST, FUNCT, JACOB, REENT,
     +   REENT1, SHOW, SOLVE, STORE, SUCCES, SUCCE1, TIME
C*****precision > double
      DOUBLE PRECISION
C*****END precision > double
C*****precision > single
C      REAL
C*****END precision > single
     +   ABOVE, BELOW, BUFFER, CHANGE, CONDIT, CSAVE, DUMMY, S, SIZE0,
     +   SIZE1, STRIAL, TDEC, TINC, TMAX, TMIN, TOLER, X, XSAVE, XTRIAL,
     +   Y, YNORM, YTRIAL

      PARAMETER (ID = 'TIMSTP:  ')
      PARAMETER (TRMAX = 10)

C     TIME STEP SIZE STATUS CODES
      PARAMETER
     +   (QBAD = 1, QDEC = 2, QGOOD = 3, QNONE = 4, QNOCH = 5, QINC = 6)

C     REPORT CODES
      PARAMETER
     +   (QNULL = 0, QDONE = 1, QDVRG = 2, QBNDS = 3, QSTDY = 4,
     +   QSMAX = 5)

      DIMENSION
     +   ABOVE(N), BELOW(N), BUFFER(N), COLUMN(7), HEADER(5, 2), S(N),
     +   STRIAL(N), X(N), XSAVE(N), XTRIAL(N), Y(N), YTRIAL(N)

C///  SAVE LOCAL VALUES DURING RETURNS FOR REVERSE COMMUNCIATION.

      SAVE

C///////////////////////////////////////////////////////////////////////
C
C     PROLOGUE.
C
C///////////////////////////////////////////////////////////////////////

C///  INITIALIZE.

C     TURN OFF ALL REVERSE COMMUNICATION FLAGS.
      DECIDE = .FALSE.
      FUNCT = .FALSE.
      JACOB = .FALSE.
      SHOW = .FALSE.
      SOLVE = .FALSE.
      STORE = .FALSE.
      TIME = .FALSE.

C     TURN OFF ALL COMPLETION STATUS FLAGS.
      ERROR = .FALSE.
      REPORT = QNULL
      SUCCES = .FALSE.

C///  IF THIS IS A RETURN CALL, THEN CONTINUE WHERE THE PROGRAM PAUSED.

      IF (REENT) THEN
         REENT = .FALSE.
         GO TO ROUTE
      END IF

C///  CHECK THE ARGUMENTS.

      ERROR = .NOT. (0 .LT. DESIRE)
      IF (ERROR) GO TO 9001

      ERROR = .NOT. (1.0 .LE. TDEC .AND. 1.0 .LE. TINC)
      IF (ERROR) GO TO 9002

      ERROR = .NOT. (0 .LT. RETIRE)
      IF (ERROR) GO TO 9003

      ERROR = .NOT. (0 .LE. STEP)
      IF (ERROR) GO TO 9004

C///  LEVEL 1 AND LEVEL 2 PRINTING.

C                     123456789_123456789_123456789_123456
C                     123456   123456   123456   123456
      HEADER(1, 1) = '                  LOG10 MAX NORM    '
      HEADER(2, 1) = '          LOG10   ---------------   '
      HEADER(3, 1) = '  TIME     TIME   STEADY   CHANGE   '
      HEADER(4, 1) = '  STEP     STEP    STATE   TO THE   '
      HEADER(5, 1) = 'NUMBER     SIZE   RESIDL   SOLUTN   '

C                     123456789_123456789_12345
C                     1234 1234 123456
      HEADER(1, 2) = 'NEWTON ALGORITHM         '
      HEADER(2, 2) = '----------------         '
      HEADER(3, 2) = 'STEPS                    '
      HEADER(4, 2) = '   JACOBIANS             '
      HEADER(5, 2) = '       CONDITION   REMARK'

      IF (0 .LT. TEXT .AND. 1 .EQ. LEVEL1) THEN
         WRITE (TEXT, 10001) ID, ((HEADER(J, K), K = 1, 2), J = 1, 5)
      ELSE IF (0 .LT. TEXT .AND. 1 .LT. LEVEL1) THEN
         WRITE (TEXT, 20001) ID
      END IF

C///  EVALUATE THE STEADY STATE RESIDUAL AT THE INITIAL SOLUTION.

      ASSIGN 1100 TO ROUTE
      GO TO 7011
1100  CONTINUE
      CALL NORM (N, YNORM, BUFFER)

C///  LEVEL 1 AND LEVEL 2 PRINTING.

      IF (0 .LT. TEXT .AND. 1 .EQ. LEVEL1) THEN
         WRITE (COLUMN(1), '(I6)') STEP
         COLUMN(2) = ' '
         CALL LOGSTR ('(F6.2)', 6, COLUMN(3), YNORM)
         WRITE (TEXT, 10002) (COLUMN(J), J = 1, 3)
      ELSE IF (0 .LT. TEXT .AND. 1 .LT. LEVEL1) THEN
         CALL LOGSTR ('(F10.2)', 10, STRING, YNORM)
         WRITE (TEXT, 20002) STRING
      END IF

C///////////////////////////////////////////////////////////////////////
C
C     TOP OF THE LOOP OVER THE TIME STEPS.
C
C///////////////////////////////////////////////////////////////////////

C///  SAVE THE LATEST SOLUTION SHOULD NEWTON FAIL; STORE IT FOR USE BY
C///  THE FUNCTION.

      CALL COPY (N, X, XSAVE)
      ASSIGN 1200 TO ROUTE
      GO TO 7021
1200  CONTINUE

C///  TOP OF THE LOOP.

      BEGIN = STEP
      IF (STEP .EQ. 0) THEN
         AGE = 0
         QMOVE = QNONE
         QSTAT = QNONE
         TRIALS = 0
      END IF

1300  CONTINUE

C///////////////////////////////////////////////////////////////////////
C
C     USE NEWTON TO SOLVE THE TIME DIFFERENCED EQUATIONS.
C
C///////////////////////////////////////////////////////////////////////

C///  LEVEL 2 PRINTING.

      IF (0 .LT. TEXT .AND. 1 .LT. LEVEL1) THEN
         CALL LOGSTR ('(F10.2)', 10, STRING, SIZE1)
         WRITE (TEXT, 20003) STRING
      END IF

C///  CALL NEWTON.

      COUNT = 0
      EXIST = (0 .LT. AGE) .AND. (BEGIN .LT. STEP)
      REENT1 = .FALSE.

2100  CONTINUE
C     SUBROUTINE NEWTON
C    +   (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
C    +   EXIST, FUNCT, JACOB, LEVEL1, LEVEL2, N, REENT, REPORT, RETIRE,
C    +   S, SHOW, SOLVE, STEPS, STMAX, STRIAL, SUCCES, X, XTRIAL, Y,
C    +   YNORM, YTRIAL)
      CALL NEWTON
     +   (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
     +   EXIST, FUNCT, JACOB, LEVEL1 - 1, LEVEL2 - 1, N, REENT1, XREPOR,
     +   TDAGE, S, SHOW, SOLVE, NUMBER, STMAX, STRIAL, SUCCE1, X,
     +   XTRIAL, Y, DUMMY, YTRIAL)
      IF (ERROR) GO TO 9006

C///  PASS MOST NEWTON REQUESTS TO THE CALLER.  TREAT JACOBIAN REQUESTS
C///  SEPARATELY TO RETAIN THE CONDITION ESTIMATE AND THE COUNT OF
C///  JACOBIANS.

      IF (REENT1) THEN
         TIME = .TRUE.
         IF (DECIDE .AND. 0 .EQ. NUMBER) THEN
C           FORCE NEWTON TO TAKE AT LEAST ONE STEP
            ACCEPT = .FALSE.
            GO TO 2100
         ELSE IF (JACOB) THEN
            GO TO 7031
         ELSE
            GO TO 7041
         END IF
      END IF

C///  DETERMINE THE CHANGE TO X AND EVALUATE THE STEADY STATE RESIDUAL.

      IF (.NOT. SUCCE1) GO TO 2400
         DO 2200 J = 1, N
            BUFFER(J) = X(J) - XSAVE(J)
2200     CONTINUE
         CALL NORM (N, CHANGE, BUFFER)

         ASSIGN 2300 TO ROUTE
         GO TO 7011
2300     CONTINUE
         CALL NORM (N, YNORM, BUFFER)
2400  CONTINUE

C///  DETERMINE THE STATUS OF THE TIME STEP.

      IF (SUCCE1) THEN
         IF (CHANGE .EQ. 0.0) THEN
            QSTAT = QNOCH
         ELSE
            QSTAT = QGOOD
            AGE = AGE + 1
            STEP = STEP + 1
         END IF
      ELSE
         QSTAT = QBAD
      END IF

C///  IF NEWTON FAILED, THEN RESTORE THE SOLUTION.

      IF (.NOT. SUCCE1) THEN
         CALL COPY (N, XSAVE, X)

C///  OTHERWISE SAVE THE LATEST SOLUTION SHOULD NEWTON FAIL; STORE IT
C///  FOR USE BY THE FUNCTION.

      ELSE
         CALL COPY (N, X, XSAVE)
         ASSIGN 2500 TO ROUTE
         GO TO 7021
      END IF
2500  CONTINUE

C///////////////////////////////////////////////////////////////////////
C
C     REPORT THE OUTCOME OF THE TIME STEP.
C
C///////////////////////////////////////////////////////////////////////

C///  LEVEL 1 PRINTING.

      IF (0 .LT. TEXT .AND. 1 .EQ. LEVEL1) THEN

      DO 3100 J = 1, 7
         COLUMN(J) = ' '
3100  CONTINUE

C     COLUMN 1: NUMBER OF THE TIME STEP
      IF (QSTAT .EQ. QGOOD) WRITE (COLUMN(1), '(I6)') STEP

C     COLUMN 2: SIZE OF THE TIME STEP
      WRITE (COLUMN(2), '(F6.2)') LOG10 (SIZE1)

C     COLUMN 3: NORM OF THE STEADY STATE RESIDUAL
      IF (QSTAT .EQ. QGOOD .OR. QSTAT .EQ. QNOCH)
     +   CALL LOGSTR ('(F6.2)', 6, COLUMN(3), YNORM)

C     COLUMN 4: NORM OF THE CHANGE TO THE SOLUTION
      IF (QSTAT .EQ. QGOOD) CALL LOGSTR ('(F6.2)', 6, COLUMN(4), CHANGE)

C     COLUMN 5: NUMBER OF NEWTON STEPS
      WRITE (COLUMN(5), '(I4)') NUMBER

C     COLUMN 6: NUMBER OF JACOBIANS
      IF (0 .LT. COUNT) WRITE (COLUMN(6), '(I4)') COUNT

C     COLUMN 7: MAXIMUM CONDITION NUMBER
      IF (0 .LT. COUNT) CALL LOGSTR ('(F6.2)', 6, COLUMN(7), CSAVE)

C     REMARK:
      IF (XREPOR .EQ. QBNDS) THEN
         REMARK = 'AT BOUNDS'
      ELSE IF (XREPOR .EQ. QDVRG) THEN
         REMARK = 'DIVERGING'
      ELSE IF (XREPOR .EQ. QSMAX) THEN
         REMARK = 'STEP LIMIT'
      ELSE IF (QSTAT .EQ. QNOCH) THEN
         REMARK = 'NO CHANGE'
      ELSE
         REMARK = ' '
      END IF
      CALL EXTENT (LENGTH, REMARK)

      WRITE (TEXT, 10002) COLUMN, REMARK (1 : LENGTH)

C///  LEVEL 2 PRINTING.

      ELSE IF (0 .LT. TEXT .AND. 1 .LT. LEVEL1) THEN
         IF (QSTAT .EQ. QGOOD .OR. QSTAT .EQ. QNOCH) THEN
            CALL LOGSTR ('(F10.2)', 10, STRING, YNORM)
            WRITE (TEXT, 20004) ID, STEP, STRING
         ELSE
            WRITE (TEXT, 20005) ID
         END IF
      END IF

C///////////////////////////////////////////////////////////////////////
C
C     DECISION TREE.
C
C///////////////////////////////////////////////////////////////////////

C///  DECIDE WHETHER TO CONTINUE.

      IF (QSTAT .EQ. QGOOD) THEN
         AGAIN = TOLER .LT. YNORM
      ELSE
         AGAIN = TRIALS .LT. TRMAX
      END IF

      IF (AGAIN) THEN

C///  CONTINUE WITH THE TIME STEP.

      IF (QSTAT .EQ. QGOOD .AND. (AGE .LT. RETIRE .OR. TINC .EQ. 1.0))
     +   THEN
         Q0 = QNONE
         QMOVE = QNONE
         TRIALS = 0

C///  CHANGE THE TIME STEP.

      ELSE
         TRIALS = TRIALS + 1

C///  INCREASE THE TIME STEP BY THE SPECIFIED AMOUNT.

      IF ((QSTAT .EQ. QGOOD .AND. AGE .EQ. RETIRE) .OR.
     +   QMOVE .EQ. QNONE .AND. QSTAT .EQ. QNOCH .AND. 1.0 .LT. TINC)
     +   THEN

         IF (0 .LT. TEXT .AND. 1 .LT. LEVEL1) WRITE (TEXT, 20006)

         AGE = 0
         QMOVE = QINC
         SIZE0 = SIZE1
         SIZE1 = MIN (TMAX, SIZE1 * TINC)

C///  REDUCE THE INCREASE BY THE SQUARE ROOT.

      ELSE IF (QMOVE .EQ. QINC .AND. QSTAT .EQ. QBAD) THEN

         IF (0 .LT. TEXT .AND. 1 .LT. LEVEL1) WRITE (TEXT, 20007)

         AGE = 0
         SIZE1 = SQRT (SIZE0 * SIZE1)

C///  REDUCE THE TIME STEP BY THE SPECIFIED AMOUNT.

      ELSE IF ((QMOVE .EQ. QDEC .OR. QMOVE .EQ. QNONE) .AND.
     +   QSTAT .EQ. QBAD .AND. 1.0 .LT. TDEC) THEN

         IF (0 .LT. TEXT .AND. 1 .LT. LEVEL1) WRITE (TEXT, 20007)

         AGE = 0
         QMOVE = QDEC
         SIZE0 = SIZE1
         SIZE1 = MAX (TMIN, SIZE1 / TDEC)

C///  REDUCE THE REDUCTION BY THE SQUARE ROOT.

      ELSE IF (QMOVE .EQ. QDEC .AND. QSTAT .EQ. QNOCH) THEN

         IF (0 .LT. TEXT .AND. 1 .LT. LEVEL1) WRITE (TEXT, 20006)

         AGE = 0
         SIZE1 = SQRT (SIZE0 * SIZE1)

C///  BOTTOM OF THE TIME STEP SELECTION TREE.

            ELSE
               ERROR = .TRUE.
               GO TO 9005
            END IF
         END IF

         IF (STEP .LT. BEGIN + DESIRE) GO TO 1300
      END IF

C///////////////////////////////////////////////////////////////////////
C
C     EPILOGUE.
C
C///////////////////////////////////////////////////////////////////////

C///  LEVEL 1 PRINTING.

      IF (0 .LT. TEXT .AND. 1 .EQ. LEVEL1) THEN

      DO 3200 J = 1, 7
         COLUMN(J) = ' '
3200  CONTINUE

C     COLUMN 1: NUMBER OF THE TIME STEP
      WRITE (COLUMN(1), '(I6)') BEGIN + DESIRE

C     COLUMN 3: NORM OF THE STEADY STATE RESIDUAL
      CALL LOGSTR ('(F6.2)', 6, COLUMN(3), TOLER)

      REMARK = 'LIMITS'
      CALL EXTENT (LENGTH, REMARK)

      WRITE (TEXT, 10002) COLUMN, REMARK (1 : LENGTH)

C///  LEVEL 2 PRINTING.

      ELSE IF (0 .LT. TEXT .AND. 1 .LT. LEVEL1) THEN
         IF (QSTAT .EQ. QGOOD) THEN
            IF (YNORM .LE. TOLER) THEN
               WRITE (TEXT, 20008)
            ELSE
               WRITE (TEXT, 20009)
            END IF
         ELSE IF (QMOVE .EQ. QNONE) THEN
            IF (QSTAT .EQ. QBAD .AND. TDEC .LE. 1.0) THEN
               WRITE (TEXT, 20010)
            ELSE IF (QSTAT .EQ. QNOCH .AND. TINC .EQ. 1.0) THEN
               WRITE (TEXT, 20011)
            END IF
         ELSE
            WRITE (TEXT, 20012)
         END IF
      END IF

C///  PRINT THE SOLUTION.

      IF (0 .LT. TEXT .AND. 1 .EQ. LEVEL2 .AND. LEVEL2 .LE. LEVEL1) THEN
         WRITE (TEXT, 20013)
         GO TO 7051
      END IF
3300  CONTINUE

C///  SET THE COMPLETION STATUS FLAGS.

      IF (STEP .EQ. BEGIN + DESIRE) THEN
         SUCCES = .TRUE.
         IF (YNORM .LE. TOLER) THEN
            REPORT = QSTDY
         ELSE
            REPORT = QDONE
         END IF
      ELSE IF (BEGIN .LT. STEP) THEN
         SUCCES = .TRUE.
         REPORT = XREPOR
      ELSE
         SUCCES = .FALSE.
         REPORT = XREPOR
      END IF

      GO TO 99999

C///////////////////////////////////////////////////////////////////////
C
C     REVERSE COMMUNICATION BLOCKS.
C
C///////////////////////////////////////////////////////////////////////

C///  EVALUATE THE STEADY STATE RESIDUAL.

7011  CONTINUE
      CALL COPY (N, X, BUFFER)
      REENT = .TRUE.
      FUNCT = .TRUE.
      TIME = .FALSE.
      GO TO 99999

C///  STORE THE LATEST SOLUTION FOR USE BY THE FUNCTION.

7021  CONTINUE
      CALL COPY (N, X, BUFFER)
      REENT = .TRUE.
      STORE = .TRUE.
      GO TO 99999

C///  PASS REQUESTS FOR JACOBIANS TO THE CALLER.

7031  CONTINUE
      REENT = .TRUE.
      ASSIGN 7032 TO ROUTE
      GO TO 99999
7032  CONTINUE
      IF (COUNT .EQ. 0) THEN
         CSAVE = CONDIT
      ELSE
         CSAVE = MAX (CONDIT, CSAVE)
      END IF
      COUNT = COUNT + 1
      GO TO 2100

C///  PASS NEWTON REQUESTS TO THE CALLER.

7041  CONTINUE
      REENT = .TRUE.
      ASSIGN 2100 TO ROUTE
      GO TO 99999

C///  PRINT THE SOLUTION.

7051  CONTINUE
      CALL COPY (N, X, BUFFER)
      REENT = .TRUE.
      SHOW = .TRUE.
      ASSIGN 3300 TO ROUTE
      GO TO 99999

C///////////////////////////////////////////////////////////////////////
C
C     INFORMATIVE MESSAGES.
C
C///////////////////////////////////////////////////////////////////////

10001 FORMAT
     +  (/1X, A9, 'PERFORM A TIME INTEGRATION.'
     +   /5(/10X, A36, A25)/)

10002 FORMAT
     +  (10X, A6, 3X, A6, 3X, A6, 3X, A6, 3X, A4, 1X, A4, 1X, A6, 3X, A)

20001 FORMAT
     +  (/1X, A9, 'BEGIN A TIME INTEGRATION.')

20002 FORMAT
     +  (/10X, A10, '  LOG10 MAX NORM STEADY STATE RESIDUAL')

20006 FORMAT
     +  (/10X, 'INCREASING THE TIME STEP SIZE.')

20007 FORMAT
     +  (/10X, 'DECREASING THE TIME STEP SIZE.')

20003 FORMAT
     +  (/10X, 'CALLING NEWTON TO ATTEMPT A TIME STEP.'
     +   //10X, A10, '  LOG10 TIME STEP SIZE')

20004 FORMAT
     +  (/1X, A9, 'NEWTON COMPLETED THE TIME STEP.'
     +   //10X, I10, '  TIME STEP NUMBER'
     +    /10X, A10, '  LOG10 MAX NORM STEADY STATE RESIDUAL')

20005 FORMAT
     +  (/1X, A9, 'NEWTON DID NOT COMPLETE THE TIME STEP.')

20008 FORMAT
     +  (/10X, 'THE STEADY STATE RESIDUAL IS SO SMALL THAT A STEADY'
     +   /10X, 'STATE SOLUTION MAY HAVE BEEN OBTAINED.')

20009 FORMAT
     +  (/10X, 'THE DESIRED NUMBER OF TIME STEPS HAS BEEN PERFORMED.')

20010 FORMAT
     +  (/10X, 'THE TIME STEP IS NOT ALLOWED TO DECREASE.')

20011 FORMAT
     +  (/10X, 'THE TIME STEP IS NOT ALLOWED TO INCREASE.')

20012 FORMAT
     +  (/10X, 'THE TIME STEP WAS NOT COMPLETED DESPITE REPEATED'
     +   /10X, 'CHANGES TO THE TIME STEP SIZE.')

20013 FORMAT
     +  (/10X, 'THE SOLUTION:')

C///////////////////////////////////////////////////////////////////////
C
C     ERROR MESSAGES.
C
C///////////////////////////////////////////////////////////////////////

9001  IF (TEXT .GT. 0) WRITE (TEXT, 99001) ID, DESIRE
      GO TO 99999

9002  IF (TEXT .GT. 0) WRITE (TEXT, 99002) ID, TDEC, TINC
      GO TO 99999

9003  IF (TEXT .GT. 0) WRITE (TEXT, 99003) ID, RETIRE
      GO TO 99999

9004  IF (TEXT .GT. 0) WRITE (TEXT, 99004) ID, STEP
      GO TO 99999

9005  IF (TEXT .GT. 0) WRITE (TEXT, 99005) ID, AGE, TRIALS, QMOVE, QSTAT
      GO TO 99999

9006  IF (TEXT .GT. 0) WRITE (TEXT, 99006) ID
      GO TO 99999

99001 FORMAT
     +  (/1X, A9, 'ERROR.  THE DESIRED NUMBER OF TIME STEPS MUST BE'
     +   /10X, 'POSITIVE.'
     +  //1X, I10, '  DESIRED NUMBER')

99002 FORMAT
     +  (/1X, A9, 'ERROR.  THE FACTORS BY WHICH THE TIME STEP MAY'
     +   /10X, 'CHANGE MUST BE GREATER THAN ONE.'
     +  //1X, E10.2, '  INCREASE FACTOR'
     +   /1X, E10.2, '  DECREASE FACTOR')

99003 FORMAT
     +  (/1X, A9, 'ERROR.  THE RETIREMENT AGE OF TIME STEP SIZES MUST'
     +   /10X, 'BE POSITIVE.'
     +  //1X, I10, '  AGE')

99004 FORMAT
     +  (/1X, A9, 'ERROR.  THE STEP NUMBER MUST BE ZERO OR POSITIVE.'
     +  //1X, I10, '  NUMBER')

99005 FORMAT
     +  (/1X, A9, 'ERROR.  THE BEHAVIOR OF THE NEWTON ALGORITHM IS NOT'
     +   /10X, 'CONSISTENT WITH THE SIZE OF THE TIME STEP.'
     +  //1X, I10, '  AGE OF THE PRESENT STEP'
     +   /1X, I10, '  TRIALS TO FIND A GOOD STEP'
     +  //1X, I10, '  QMOVE'
     +   /1X, I10, '  QSTAT')

99006 FORMAT
     +  (/1X, A9, 'ERROR.  NEWTON ERRS.')

C///  EXIT.

99999 CONTINUE
      RETURN
      END
      SUBROUTINE TIMWOR (VALUE, WORD)

      IMPLICIT COMPLEX (A - Z)
C     CHARACTER DIGIT*1, WORD*9
C     THE CRAY/CTSS FORTLIB I/O ROUTINES REQUIRE THAT CHARACTER STRINGS
C     USED AS FORMATS OR INTERNAL FILES HAVE WORD BOUNDARIES.  THUS, THE
C     STRINGS' LENGTHS MUST BE DIVISIBLE BY 8.
      CHARACTER DIGIT*1, WORD*16
C*****precision > double
      DOUBLE PRECISION SECOND, SIXTY, VALUE, ZERO
C*****END precision > double
      INTEGER HOUR, J, MINUTE
      LOGICAL FOUND
C*****precision > single
C      REAL SECOND, SIXTY, VALUE, ZERO
C*****END precision > single

C     PRECISION-INDEPENDENT CONSTANT
      ZERO = 0.0
      SIXTY = 60.0

      WORD = ' '
      IF (VALUE .GT. ZERO) THEN
         SECOND = MOD (VALUE, SIXTY)
         MINUTE = MOD (INT (VALUE) / 60, 60)
         HOUR = INT (VALUE) / 3600

         IF (HOUR .LE. 9) THEN
            WRITE (WORD, 10001, ERR = 0200) HOUR, MINUTE, SECOND
         ELSE
C     THE LINE BELOW WAS REPLACED WITH A FUNCTIONALLY EQUAIVALENT LINE
C     BECAUSE CRAY/CTSS/FORTLIB LACKS THE INTRINSIC FUNCTION IDNINT.
C           WRITE (WORD, 10002, ERR = 0200) HOUR, MINUTE, INT(SECOND)
            WRITE (WORD, 10002, ERR = 0200)
     +         HOUR, MINUTE, NINT (REAL (SECOND))
         END IF

         FOUND = .FALSE.
         DO 0100 J = 1, 9
            DIGIT = WORD(J : J)
            FOUND = FOUND .OR. (DIGIT .NE. ' ' .AND. DIGIT .NE. '0'
     +         .AND. DIGIT .NE. ':')
            IF (FOUND) THEN
               IF (DIGIT .EQ. ' ') WORD(J : J) = '0'
            ELSE
               WORD(J : J) = ' '
            END IF
0100     CONTINUE
      END IF
0200  CONTINUE

10001 FORMAT (I1, ':', I2, ':', F4.1)
10002 FORMAT (I3, ':', I2, ':', I2)

      RETURN
      END
      SUBROUTINE TWOPNT
     +  (ERROR, TEXT, LEVEL,
     +   ISIZE, IWORK, RSIZE, RWORK,
     +   ABOVE, ACTIVE, ADAPT, BELOW, BINARY, BUFFER, COMPS, CONDIT,
     +   FUNCT, JACOB, MARK, MESH, PASS, PASSES, PLIMIT, PMAX, POINTS,
     +   REENT, REPORT, SAVE, SCALRS, SHOW, SOLVE, SSABS, SSAGE, SSREL,
     +   STEPS0, STEPS1, STEPS2, STORE, SUCCES, TDABS, TDAGE, TDEC,
     +   TDREL, TIME, TINC, TMAX, TMIN, TOLER0, TOLER1, TOLER2, TSTEP0,
     +   TSTEP1, UPDATE, X)

      IMPLICIT COMPLEX (A - P, R - Z), INTEGER (Q)
      CHARACTER
     +   BANNER*60, COLUMN*16, HEADER*60, ID*9, LINE*60, REMARK*60,
     +   WORD*16
      EXTERNAL
     +   COPY, CPUTIM, ELAPSE, EXTENT, GRAB, LOGSTR, NEWTON, NORM,
     +   REFINE, TIMSTP, TIMWOR
      INTEGER
     +   BINARY, COMPS, DESIRE, EVENT, GMAX, GRID, ISIZE, IWORK, J,
     +   JACOBS, K, LABEL, LAST, LENGTH, LEVEL, PASS, PASSES, PLIMIT,
     +   PMAX, POINTS, REPORT, RETURN, ROUTE, RSIZE, SCALRS, SIZE,
     +   SSAGE, STEP, STEPS, STEPS0, STEPS1, STEPS2, STMAX, TDAGE, TEXT,
     +   XREPOR
      LOGICAL
     +   ACCEPT, ACTIVE, ADAPT, DECIDE, ERROR, EXIST, FIRST, FOUND,
     +   FUNCT, JACOB, MARK, REENT, SATISF, SAVE, SHOW, SOLVE, STORE,
     +   TIME, UPDATE, XREENT, SUCCES
C*****precision > double
      DOUBLE PRECISION
C*****END precision > double
C*****precision > single
C      REAL
C*****END precision > single
     +   ABOVE, BELOW, BUFFER, CENT, CONDIT, DETAIL, MAXCON, MESH,
     +   RWORK, SSABS, SSREL, TDABS, TDEC, TDREL, TEMP, TIMER, TINC,
     +   TMAX, TMIN, TOLER, TOLER0, TOLER1, TOLER2, TOTAL, TSTEP0,
     +   TSTEP1, X, YNORM, ZERO

      PARAMETER (ID = 'TWOPNT:  ')
      PARAMETER (GMAX = 100)
      PARAMETER (TOLER = 1.0E-10)

C     REPORT CODES
      PARAMETER (QNULL = 0)

C     REPORT CODES FROM NEWTON AND TIMSTP
      PARAMETER (QDONE = 1, QDVRG = 2, QBNDS = 3, QSTDY = 4, QSMAX = 5)

C     VALUES OF REPORT UPON EXIT WITHOUT SUCCESS AND WITHOUT ERROR.
      PARAMETER (QNSLVE = 1, QNSPAC = 2)

C     LOCATION OF DATA IN ARRAYS DETAIL, EVENT, TIMER, AND TOTAL.  THE
C     LOCATIONS ARE CHOSEN TO SIMPLIFY WRITE STATEMENTS.  DETAIL USES
C     ONLY 1 : 8, EVENT USES ONLY 5 : 8, TIMER USES 1 : 9, AND TOTAL
C     USES ONLY 2 : 9.  IN ADDITION, 2, 3, 4, 10, AND 11 ARE USED AS
C     MNEMONIC VALUES FOR QTASK.
      PARAMETER
     +   (QGRID  =  1,
     +    QTIMST =  2,
     +    QNEWTO =  3,
     +    QREFIN =  4,
     +    QFUNCT =  5,
     +    QJACOB =  6,
     +    QSOLVE =  7,
     +    QOTHER =  8,
     +    QTOTAL =  9,
     +    QENTRY = 10,
     +    QEXIT  = 11)

      DIMENSION
     +   ABOVE(SCALRS + COMPS * PMAX), ACTIVE(COMPS), BANNER(4),
     +   BELOW(SCALRS + COMPS * PMAX), BUFFER(SCALRS + COMPS * PMAX),
     +   COLUMN(6), DETAIL(GMAX, QTOTAL), EVENT(GMAX, QTOTAL),
     +   HEADER(3, 3), IWORK(ISIZE), LEVEL(2), MARK(PMAX), MESH(PMAX),
     +   RWORK(RSIZE), SIZE(GMAX), TIMER(QTOTAL), TOTAL(QTOTAL),
     +   X(SCALRS + COMPS * PMAX)

C///  SAVE LOCAL VALUES DURING RETURNS FOR REVERSE COMMUNCIATION.

      SAVE

C///////////////////////////////////////////////////////////////////////
C
C     RETURN FROM REVERSE COMMUNICATION.
C
C///////////////////////////////////////////////////////////////////////

C///  TURN OFF ALL REVERSE COMMUNICATION FLAGS.

      FUNCT = .FALSE.
      JACOB = .FALSE.
      SAVE = .FALSE.
      SHOW = .FALSE.
      SOLVE = .FALSE.
      STORE = .FALSE.
      TIME = .FALSE.
      UPDATE = .FALSE.

C///  TURN OFF THE ERROR FLAG.

      ERROR = .FALSE.

C///  CONTINUE WHERE THE PROGRAM PAUSED.

      IF (REENT) THEN
         REENT = .FALSE.
         GO TO ROUTE
      END IF

C///  TURN OFF THE COMPLETION STATUS FLAGS.

      REPORT = QNULL
      SUCCES = .FALSE.

C///////////////////////////////////////////////////////////////////////
C
C     ENTRY BLOCK.  INITIALIZE A NEW PROBLEM.
C
C///////////////////////////////////////////////////////////////////////

C///  CHECK THE ARGUMENTS.

      ERROR = LEVEL(1) .LT. LEVEL(2)
      IF (ERROR) GO TO 9001

      ERROR = .NOT. (0 .LT. PASSES)
      IF (ERROR) GO TO 9002

      ERROR = .NOT. (0 .LE. SCALRS .AND. 0 .LE. COMPS
     +   .AND. 0 .LE. POINTS .AND. POINTS .LE. PMAX
     +   .AND. ((0 .LT. COMPS) .EQV. (0 .LT. POINTS)))
      IF (ERROR) GO TO 9003

      ERROR = .NOT.
     +   (1 .LE. COMPS .AND. 1 .LE. PMAX .AND.
     +   1 .LE. SCALRS + COMPS * PMAX)
      IF (ERROR) GO TO 9004

C///  ONE-TIME INITIALIZATION OF CONSTANTS AND VARIABLES.

C     DECISION BLOCK PARAMETERS.
      QTASK = QENTRY

C     PRECISION-INDEPENDENT CONSTANTS.
      ZERO = 0
      CENT = 100

C     STATISTICS POINTER AND ARRAYS.
      GRID = 1
      DO 1100 K = 1, QTOTAL
         TOTAL(K) = ZERO
         DO 1100 J = 1, GMAX
            DETAIL(J, K) = ZERO
            EVENT(J, K) = 0
1100  CONTINUE

C     MARKERS FOR NEW POINTS.
      DO 1200 K = 1, POINTS
         MARK(K) = .FALSE.
1200  CONTINUE

C     TIME STEP NUMBER AND SIZE.
      STEP = 0
      TSTEP1 = TSTEP0

C///  PARTITION THE INTEGER WORK SPACE.

C     SUBROUTINE GRAB (ERROR, TEXT, LAST, FIRST, NUMBER)

      LAST = 0

C     VARY(PMAX)
      CALL GRAB (ERROR, TEXT, LAST, QVARY, PMAX)
      IF (ERROR) GO TO 9005

C     VARY1(PMAX)
      CALL GRAB (ERROR, TEXT, LAST, QVARY1, PMAX)
      IF (ERROR) GO TO 9005

C     VARY2(PMAX)
      CALL GRAB (ERROR, TEXT, LAST, QVARY2, PMAX)
      IF (ERROR) GO TO 9005

      ERROR = .NOT. (LAST .LE. ISIZE)
      IF (ERROR) GO TO 9006

C///  PARTITION THE REAL WORK SPACE.

C     SUBROUTINE GRAB (ERROR, TEXT, LAST, FIRST, NUMBER)

      LAST = 0

C     RATIO1(PMAX)
      CALL GRAB (ERROR, TEXT, LAST, QRAT1, PMAX)
      IF (ERROR) GO TO 9005

C     RATIO2(PMAX)
      CALL GRAB (ERROR, TEXT, LAST, QRAT2, PMAX)
      IF (ERROR) GO TO 9005

C     S(SCALRS + COMPS * PMAX)
      CALL GRAB (ERROR, TEXT, LAST, QS, SCALRS + COMPS * PMAX)
      IF (ERROR) GO TO 9005

C     STRIAL(SCALRS + COMPS * PMAX)
      CALL GRAB (ERROR, TEXT, LAST, QSTRL, SCALRS + COMPS * PMAX)
      IF (ERROR) GO TO 9005

C     XSAVE(SCALRS + COMPS * PMAX)
      CALL GRAB (ERROR, TEXT, LAST, QXSAV, SCALRS + COMPS * PMAX)
      IF (ERROR) GO TO 9005

C     XTEMP(SCALRS + COMPS * PMAX)
      CALL GRAB (ERROR, TEXT, LAST, QXTMP, SCALRS + COMPS * PMAX)
      IF (ERROR) GO TO 9005

C     XTRIAL(SCALRS + COMPS * PMAX)
      CALL GRAB (ERROR, TEXT, LAST, QXTRL, SCALRS + COMPS * PMAX)
      IF (ERROR) GO TO 9005

C     Y(SCALRS + COMPS * PMAX)
      CALL GRAB (ERROR, TEXT, LAST, QY, SCALRS + COMPS * PMAX)
      IF (ERROR) GO TO 9005

C     YTRIAL(SCALRS + COMPS * PMAX)
      CALL GRAB (ERROR, TEXT, LAST, QYTRL, SCALRS + COMPS * PMAX)
      IF (ERROR) GO TO 9005

      ERROR = .NOT. (LAST .LE. RSIZE)
      IF (ERROR) GO TO 9007

C///  INITIALIZE THE TOTAL TIME STATISTIC.

      CALL CPUTIM (TIMER(QTOTAL))

C///  INITIALIZE THE STATISTICS FOR THE FIRST MESH.

      IF (GRID .LE. GMAX) SIZE(GRID) = POINTS
      IF (GRID .LE. GMAX) CALL CPUTIM (TIMER(QGRID))

C///  SAVE THE INITIAL SOLUTION.

      ASSIGN 1300 TO RETURN
      GO TO 9911
1300  CONTINUE

C///  PRINT THE BANNER AT ALL PRINT LEVELS.

      BANNER(1) = 'TWO-POINT BOUNDARY VALUE PROBLEM SOLVER,'
      BANNER(2) =
     +   'VERSION 2.22B OF MARCH 1991 BY DR. JOSEPH F. GRCAR.'
C*****precision > double
      BANNER(3) = 'DOUBLE PRECISION'
C*****END precision > double
C*****precision > single
C      BANNER(3) = 'SINGLE PRECISION'
C*****END precision > single
      IF (ADAPT) THEN
         BANNER(4) = 'PERMIT MESH REFINEMENT'
      ELSE
         BANNER(4) = 'DO NOT REFINE THE MESH'
      END IF

      IF (LEVEL(1) .GE. 1) THEN
         IF (0 .LT. TEXT) WRITE (TEXT, 10001) ID, BANNER
      END IF

C///  PRINT THE INITIAL GUESS AT PRINT LEVELS 11, 21, AND 22.

      IF (LEVEL(2) .GE. 1) THEN
         IF (0 .LT. TEXT) WRITE (TEXT, 10002) ID
         ASSIGN 1400 TO RETURN
         GO TO 9921
      END IF
1400  CONTINUE

C///  PRINT LEVEL 10 OR 11.

C                     123456789_123456789_123456789_123456
C                     123456789   123456789   123456789
      HEADER(1, 1) = '                LOG10   LOG10 MAX   '
      HEADER(2, 1) = '             MAX NORM   CONDITION   '
      HEADER(3, 1) = ' ACTIVITY    RESIDUAL      NUMBER   '

C                     123456789_123456789_123456
C                     123456789_1234   123456789
      HEADER(1, 2) = 'POINTS                    '
      HEADER(2, 2) = '    STEPS                 '
      HEADER(3, 2) = '     JACOBIANS   REMARK   '

      IF (LEVEL(1) .NE. 1) GO TO 1600
         IF (0 .LT. TEXT) WRITE (TEXT, 10003)
     +      ID, ((HEADER(J, K), K = 1, 2), J = 1, 3)
         ASSIGN 1500 TO LABEL
         GO TO 7100
1500     CONTINUE
         CALL EXTENT (LENGTH, REMARK)
         IF (0 .LT. TEXT) WRITE (TEXT, 10017)
     +      COLUMN, REMARK (1 : LENGTH)
1600  CONTINUE

C///////////////////////////////////////////////////////////////////////
C
C     DECISION BLOCK.  THE PREVIOUS TASK DETERMINES THE NEXT.
C
C///////////////////////////////////////////////////////////////////////

2100  CONTINUE

C///  ENTRY WAS THE PREVIOUS TASK.

      IF (QTASK .EQ. QENTRY) THEN
         PASS = 1
         IF (0 .LT. ABS (STEPS0)) THEN
            QTASK = QTIMST
            DESIRE = ABS (STEPS0)
         ELSE
            QTASK = QNEWTO
         END IF

C///  NEWTON WAS THE PREVIOUS TASK.

      ELSE IF (QTASK .EQ. QNEWTO) THEN
         IF (FOUND) THEN
            IF (PASS .LT. PASSES) THEN
               PASS = PASS + 1
               QTASK = QNEWTO
            ELSE
               PASS = 1
               IF (ADAPT) THEN
                  QTASK = QREFIN
               ELSE
                  QTASK = QEXIT
                  SUCCES = .TRUE.
               END IF
            END IF
         ELSE
            IF (STEPS1 .GT. 0) THEN
               QTASK = QTIMST
               DESIRE = STEPS1
            ELSE
               QTASK = QEXIT
               REPORT = QNSLVE
               SUCCES = .FALSE.
            END IF
         END IF

C///  REFINE WAS THE PREVIOUS TASK.

      ELSE IF (QTASK .EQ. QREFIN) THEN
         IF (FOUND) THEN
            PASS = 1
            STEP = 0
C           RESET THE TIME STEP.
            TSTEP1 = TSTEP0
            QTASK = QNEWTO
         ELSE
            QTASK = QEXIT
            IF (.NOT. SATISF) REPORT = QNSPAC
            SUCCES = SATISF
         END IF

C///  TIMSTP WAS THE PREVIOUS TASK.

      ELSE IF (QTASK .EQ. QTIMST) THEN
         IF (FOUND) THEN
            IF (STEPS0 .LT. 0) THEN
               QTASK = QEXIT
               SUCCES = .TRUE.
            ELSE
               QTASK = QNEWTO
            END IF
         ELSE
            QTASK = QEXIT
            REPORT = QNSLVE
            SUCCES = .FALSE.
         END IF
      END IF

C///  BRANCH TO THE NEXT TASK.

      IF (QTASK .EQ. QEXIT) GO TO 3100
      IF (QTASK .EQ. QNEWTO) GO TO 4100
      IF (QTASK .EQ. QREFIN) GO TO 5100
      IF (QTASK .EQ. QTIMST) GO TO 6100
      ERROR = .TRUE.
      GO TO 9008

C///////////////////////////////////////////////////////////////////////
C
C     EXIT BLOCK.
C
C///////////////////////////////////////////////////////////////////////

3100  CONTINUE

C///  COMPLETE STATISTICS FOR THE LAST MESH.

      CALL ELAPSE (TIMER(QGRID))
      IF (GRID .LE. GMAX) DETAIL(GRID, QGRID) = TIMER(QGRID)

C///  PRINT LEVEL 11 OR 21.

      IF (LEVEL(2) .EQ. 1) THEN
         IF (0 .LT. TEXT) WRITE (TEXT, 10005) ID
         ASSIGN 3200 TO RETURN
         GO TO 9921
      END IF
3200  CONTINUE

C///  COMPLETE THE TOTAL TIME STATISTIC.

      CALL ELAPSE (TIMER(QTOTAL))
      TOTAL(QTOTAL) = TIMER(QTOTAL)

C///  REPORT THE TOTAL TIME.

      IF (LEVEL(1) .GE. 1) THEN
         IF (TOTAL(QTOTAL) .GE. 36000.0) THEN
            LINE = 'HOURS : MINUTES : SECONDS'
         ELSE IF (TOTAL(QTOTAL) .GE. 3600.0) THEN
            LINE = 'HOURS : MINUTES : SECONDS . FRACTION'
         ELSE IF (TOTAL(QTOTAL) .GE. 60.0) THEN
            LINE = 'MINUTES : SECONDS . FRACTION'
         ELSE
            LINE = 'SECONDS . FRACTION'
         END IF
         IF (TOTAL(QTOTAL) .GT. ZERO) THEN
            CALL TIMWOR (TOTAL(QTOTAL), WORD)
            IF (0 .LT. TEXT) WRITE (TEXT, 10006) ID, WORD, LINE
         END IF
      END IF

C///  REPORT THE PERCENT OF TOTAL CPU TIME.

C                     123456789_123456
C                     123456  123456
      HEADER(1, 1) = '                '
      HEADER(2, 1) = '  GRID    GRID  '
      HEADER(3, 1) = 'POINTS   TOTAL  '

C                     123456789_123456789_12
C                     123456 123456 123456
      HEADER(1, 2) = 'ACTIVITY              '
      HEADER(2, 2) = '--------------------  '
      HEADER(3, 2) = 'TIMSTP NEWTON REFINE  '

C                     123456789_123456789_1234567
C                     123456 123456 123456 123456
      HEADER(1, 3) = 'EVENT                      '
      HEADER(2, 3) = '---------------------------'
      HEADER(3, 3) = 'FUNCTN JACOBN  SOLVE  OTHER'

      IF (LEVEL(1) .GE. 1) THEN
         IF (TOTAL(QTOTAL) .GT. ZERO) THEN
            IF (0 .LT. TEXT) WRITE (TEXT, 10007)
     +         ID, ((HEADER(J, K), K = 1, 3), J = 1, 3),
     +         (SIZE(J), (CENT * DETAIL(J, K) / TOTAL(QTOTAL),
     +         K = 1, 8), J = 1, GRID)
            IF (GRID .GT. 1) THEN
               IF (0 .LT. TEXT) WRITE (TEXT, 10008)
     +            (CENT * TOTAL(K) / TOTAL(QTOTAL), K = 2, 8)
            END IF
            IF (GRID .GT. GMAX) THEN
               IF (0 .LT. TEXT) WRITE (TEXT, 10009)
            END IF
         END IF
      END IF

C///  REPORT THE AVERAGE CPU TIME.

C                     123456789
C                     123456
      HEADER(1, 1) = '         '
      HEADER(2, 1) = '  GRID   '
      HEADER(3, 1) = 'POINTS   '

C                     123456789_123456789_12345678
C                     1234567  1234567  1234567
      HEADER(1, 2) = 'AVERAGE SECONDS             '
      HEADER(2, 2) = '-------------------------   '
      HEADER(3, 2) = ' FUNCTN   JACOBN    SOLVE   '

C                     123456789_123456789_12345
C                     1234567  1234567  1234567
      HEADER(1, 3) = 'NUMBER OF EVENTS         '
      HEADER(2, 3) = '-------------------------'
      HEADER(3, 3) = ' FUNCTN   JACOBN    SOLVE'

      IF (LEVEL(1) .GE. 1) THEN
         IF (TOTAL(QTOTAL) .GT. ZERO) THEN
            IF (0 .LT. TEXT) WRITE (TEXT, 10010)
     +         ID, ((HEADER(J, K), K = 1, 3), J = 1, 3),
     +         (SIZE(J), (DETAIL(J, K) / EVENT(J, K), K = 5, 7),
     +         (EVENT(J, K), K = 5, 7), J = 1, GRID)
            IF (GRID .GT. GMAX) THEN
               IF (0 .LT. TEXT) WRITE (TEXT, 10011)
            END IF
         END IF
      END IF

C///  REPORT THE COMPLETION STATUS.

      IF (LEVEL(1) .GE. 1) THEN
         IF (SUCCES) THEN
            IF (0 .LT. TEXT) WRITE (TEXT, 10012) ID
         ELSE IF (REPORT .EQ. QNSLVE) THEN
            IF (GRID .EQ. 1) THEN
               IF (0 .LT. TEXT) WRITE (TEXT, 10013) ID
            ELSE
               IF (0 .LT. TEXT) WRITE (TEXT, 10014) ID
            END IF
         ELSE IF (REPORT .EQ. QNSPAC) THEN
            IF (0 .LT. TEXT) WRITE (TEXT, 10015) ID
         ELSE
            ERROR = .TRUE.
            GO TO 9009
         END IF
      END IF

C///  BOTTOM OF THE EXIT BLOCK.

      GO TO 99999

C///////////////////////////////////////////////////////////////////////
C
C     NEWTON BLOCK.
C
C///////////////////////////////////////////////////////////////////////

4100  CONTINUE

C///  INITIALIZE STATISTICS ON ENTRY TO THE NEWTON BLOCK.

      CALL CPUTIM (TIMER(QNEWTO))
      FIRST = .TRUE.
      JACOBS = 0
      MAXCON = ZERO

C///  PRINT LEVEL 20, 21, OR 22 ON ENTRY TO THE NEWTON BLOCK.

      IF (LEVEL(1) .GE. 2) THEN
         IF (0 .LT. TEXT) WRITE (TEXT, 10016) ID
      END IF

C///  PREPARE TO CALL NEWTON.

C     SAVE THE SOLUTION SHOULD NEWTON FAIL.
      CALL COPY (SCALRS + COMPS * POINTS, X, RWORK(QXSAV))

      EXIST = .FALSE.

C///  CALL NEWTON.

      STMAX = 0
      XREENT = .FALSE.
4200  CONTINUE
C     SUBROUTINE NEWTON
C    +   (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
C    +   EXIST, FUNCT, JACOB, LEVEL1, LEVEL2, N, REENT, REPORT, RETIRE,
C    +   S, SHOW, SOLVE, STEPS, STMAX, STRIAL, SUCCES, X, XTRIAL, Y,
C    +   YNORM, YTRIAL)
      CALL NEWTON
     +   (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
     +   EXIST, FUNCT, JACOB, LEVEL(1) - 1, LEVEL(2) - 1,
     +   SCALRS + COMPS * POINTS, XREENT, XREPOR, SSAGE, RWORK(QS),
     +   SHOW, SOLVE, STEPS, STMAX, RWORK(QSTRL), FOUND, X,
     +   RWORK(QXTRL), RWORK(QY), YNORM, RWORK(QYTRL))
      IF (ERROR) GO TO 9010

C///  SERVICE REQUESTS FROM NEWTON: DECIDE WHETHER THE SOLUTION IS
C///  ACCEPTABLE.

      IF (XREENT) THEN
         IF (DECIDE) THEN
            ACCEPT = .TRUE.
            DO 4300 J = 1, SCALRS + COMPS * POINTS
               ACCEPT = ACCEPT .AND.
     +            ABS (RWORK(QS - 1 + J)) .LE.
     +            MAX (SSABS, SSREL * ABS (X(J)))
4300        CONTINUE
            GO TO 4200

C///  PASS OTHER REQUESTS FROM NEWTON TO THE CALLER.

         ELSE
            ASSIGN 4200 TO RETURN
            GO TO 9931
         END IF
      END IF

C///  REACT TO THE COMPLETION OF NEWTON.

      IF (FOUND) THEN
C        SAVE THE LATEST SOLUTION.
         ASSIGN 4400 TO RETURN
         GO TO 9911
      ELSE
C        RESTORE THE SOLUTION.
         CALL COPY (SCALRS + COMPS * POINTS, RWORK(QXSAV), X)
      END IF
4400  CONTINUE

C///  COMPLETE STATISTICS FOR THE NEWTON BLOCK.

      CALL ELAPSE (TIMER(QNEWTO))
      TOTAL(QNEWTO) = TOTAL(QNEWTO) + TIMER(QNEWTO)
      IF (GRID .LE. GMAX)
     +   DETAIL(GRID, QNEWTO) = DETAIL(GRID, QNEWTO) + TIMER(QNEWTO)

C///  PRINT LEVEL 10 OR 11 ON EXIT FROM THE NEWTON BLOCK.

      IF (LEVEL(1) .NE. 1) GO TO 4600
         ASSIGN 4500 TO LABEL
         GO TO 7100
4500     CONTINUE
         CALL EXTENT (LENGTH, REMARK)
         IF (0 .LT. TEXT) WRITE (TEXT, 10017)
     +      COLUMN, REMARK (1 : LENGTH)
4600  CONTINUE

C///  PRINT LEVEL 20, 21, OR 22 ON EXIT FROM THE NEWTON BLOCK.

      IF (LEVEL(1) .GE. 2) THEN
         IF (FOUND) THEN
            IF (0 .LT. TEXT) WRITE (TEXT, 10018) ID
         ELSE
            IF (0 .LT. TEXT) WRITE (TEXT, 10019) ID
         END IF
      END IF

C///  BOTTOM OF THE NEWTON BLOCK.

      GO TO 2100

C///////////////////////////////////////////////////////////////////////
C
C     REFINE BLOCK.
C
C///////////////////////////////////////////////////////////////////////

5100  CONTINUE

C///  INITIALIZE STATISTICS ON ENTRY TO THE REFINE BLOCK.

      CALL CPUTIM (TIMER(QREFIN))

C///  PRINT LEVEL 20, 21, OR 22 ON ENTRY TO THE REFINE BLOCK.

      IF (LEVEL(1) .GE. 2) THEN
         IF (0 .LT. TEXT) WRITE (TEXT, 10020) ID
      END IF

C///  CALL REFINE.

      XREENT = .FALSE.
5200  CONTINUE
C     SUBROUTINE REFINE
C    +   (ERROR, TEXT, ABOVE, ACTIVE, BELOW, BUFFER, COMPS, LEVEL1,
C    +   LEVEL2, MARK, MESH, N, NEWMSH, PLIMIT, PMAX, POINTS, RATIO1,
C    +   RATIO2, REENT, SHOW, SUCCES, TOLER0, TOLER1, TOLER2, UPDATE,
C    +   VARY, VARY1, VARY2, X)
      CALL REFINE
     +   (ERROR, TEXT, ABOVE(SCALRS + 1), ACTIVE, BELOW(SCALRS + 1),
     +   BUFFER(SCALRS + 1), COMPS, LEVEL(1) - 1, LEVEL(2) - 1, MARK,
     +   MESH, FOUND, PLIMIT, PMAX, POINTS, RWORK(QRAT1), RWORK(QRAT2),
     +   XREENT, SHOW, SATISF, TOLER0, TOLER1, TOLER2, UPDATE,
     +   IWORK(QVARY), IWORK(QVARY1), IWORK(QVARY2), X(SCALRS + 1))
      IF (ERROR) GO TO 9011

C///  SERVICE REQUESTS FROM REFINE: PASS REQUESTS TO THE CALLER.

      IF (XREENT) THEN
         IF (0 .LT. SCALRS) CALL COPY (SCALRS, X, BUFFER)
         ASSIGN 5200 TO RETURN
         GO TO 9931
      END IF

C///  REACT TO THE COMPLETION OF REFINE.

      IF (.NOT. FOUND) GO TO 5400

C        COMPLETE STATISTICS FOR THE OLD MESH.
         CALL ELAPSE (TIMER(QGRID))
         IF (GRID .LE. GMAX) DETAIL(GRID, QGRID) = TIMER(QGRID)

C        INITIALIZE STATISTICS FOR THE NEW MESH.
         GRID = GRID + 1
         IF (GRID .LE. GMAX) THEN
            CALL CPUTIM (TIMER(QGRID))
            SIZE(GRID) = POINTS
         END IF

C        SAVE THE LATEST SOLUTION.
         ASSIGN 5300 TO RETURN
         GO TO 9911
5300     CONTINUE

5400  CONTINUE

C///  COMPLETE STATISTICS FOR THE REFINE BLOCK.

      CALL ELAPSE (TIMER(QREFIN))
      TOTAL(QREFIN) = TOTAL(QREFIN) + TIMER(QREFIN)
      IF (GRID .LE. GMAX)
     +   DETAIL(GRID, QREFIN) = DETAIL(GRID, QREFIN) + TIMER(QREFIN)

C///  PRINT LEVEL 10 OR 11 ON EXIT FROM THE REFINE BLOCK.

      IF (LEVEL(1) .NE. 1) GO TO 5600
         ASSIGN 5500 TO LABEL
         GO TO 7100
5500     CONTINUE
         CALL EXTENT (LENGTH, REMARK)
         IF (0 .LT. TEXT) WRITE (TEXT, 10004)
     +      COLUMN, REMARK (1 : LENGTH)
5600  CONTINUE

C///  PRINT LEVEL 20, 21, OR 22 ON EXIT FROM THE REFINE BLOCK.

      IF (LEVEL(1) .GE. 2) THEN
         IF (FOUND) THEN
            IF (0 .LT. TEXT) WRITE (TEXT, 10021) ID
         ELSE
            IF (0 .LT. TEXT) WRITE (TEXT, 10022) ID
         END IF
      END IF
5700  CONTINUE

C///  BOTTOM OF THE REFINE BLOCK.

      GO TO 2100

C///////////////////////////////////////////////////////////////////////
C
C     TIMSTP BLOCK.
C
C///////////////////////////////////////////////////////////////////////

6100  CONTINUE

C///  INITIALIZE STATISTICS ON ENTRY TO THE TIMSTP BLOCK.

      CALL CPUTIM (TIMER(QTIMST))
      FIRST = .TRUE.
      JACOBS = 0
      MAXCON = ZERO
      STEPS = STEP

C///  PRINT LEVEL 20, 21, OR 22 ON ENTRY TO THE TIMSTP BLOCK.

      IF (LEVEL(1) .GE. 2) THEN
         IF (0 .LT. TEXT) WRITE (TEXT, 10023) ID
      END IF

C///  CALL TIMSTP.

      STMAX = 100
      XREENT = .FALSE.
6200  CONTINUE
C     SUBROUTINE TIMSTP
C    +   (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
C    +   DESIRE, FUNCT, JACOB, LEVEL1, LEVEL2, N, REENT, REPORT,
C    +   RETIRE, S, SHOW, SIZE1, SOLVE, STEP, STMAX, STORE, STRIAL,
C    +   SUCCES, TDAGE, TDEC, TIME, TINC, TMAX, TMIN, TOLER, X, XSAVE,
C    +   XTRIAL, Y, YNORM, YTRIAL)
      CALL TIMSTP
     +   (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
     +   DESIRE, FUNCT, JACOB, LEVEL(1) - 1, LEVEL(2) - 1,
     +   SCALRS + COMPS * POINTS, XREENT, XREPOR, STEPS2, RWORK(QS),
     +   SHOW, TSTEP1, SOLVE, STEP, STMAX, STORE, RWORK(QSTRL), FOUND,
     +   TDAGE, TDEC, TIME, TINC, TMAX, TMIN, TOLER, X, RWORK(QXTMP),
     +   RWORK(QXTRL), RWORK(QY), YNORM, RWORK(QYTRL))
      IF (ERROR) GO TO 9012

C///  SERVICE REQUESTS FROM TIMSTP: DECIDE WHETHER THE SOLUTION IS
C///  ACCEPTABLE.

      IF (XREENT) THEN
         IF (DECIDE) THEN
            ACCEPT = .TRUE.
            DO 6300 J = 1, SCALRS + COMPS * POINTS
               ACCEPT = ACCEPT .AND.
     +            ABS (RWORK(QS - 1 + J)) .LE.
     +            MAX (TDABS, TDREL * ABS (X(J)))
6300        CONTINUE
            GO TO 6200

C///  PASS OTHER REQUESTS FROM TIMSTP TO THE CALLER.

         ELSE
            ASSIGN 6200 TO RETURN
            GO TO 9931
         END IF
      END IF

C///  REACT TO THE COMPLETION OF TIMSTP.

      IF (FOUND) THEN
C        SAVE THE LATEST SOLUTION.
         ASSIGN 6400 TO RETURN
         GO TO 9911
      END IF
6400  CONTINUE

C///  COMPLETE STATISTICS FOR THE TIMSTP BLOCK.

      CALL ELAPSE (TIMER(QTIMST))
      TOTAL(QTIMST) = TOTAL(QTIMST) + TIMER(QTIMST)
      IF (GRID .LE. GMAX)
     +   DETAIL(GRID, QTIMST) = DETAIL(GRID, QTIMST) + TIMER(QTIMST)
      STEPS = STEP - STEPS

C///  PRINT LEVEL 10 OR 11 ON EXIT FROM THE TIMSTP BLOCK.

      IF (LEVEL(1) .NE. 1) GO TO 6600
         ASSIGN 6500 TO LABEL
         GO TO 7100
6500     CONTINUE
         CALL EXTENT (LENGTH, REMARK)
         IF (0 .LT. TEXT) WRITE (TEXT, 10017)
     +      COLUMN, REMARK (1 : LENGTH)
6600  CONTINUE

C///  PRINT LEVEL 20, 21, OR 22 ON EXIT FROM THE TIMSTP BLOCK.

      IF (LEVEL(1) .GE. 2) THEN
         IF (FOUND) THEN
            IF (0 .LT. TEXT) WRITE (TEXT, 10024) ID
         ELSE
            IF (0 .LT. TEXT) WRITE (TEXT, 10025) ID
         END IF
      END IF

C///  BOTTOM OF THE TIMSTP BLOCK.

      GO TO 2100

C///////////////////////////////////////////////////////////////////////
C
C     BLOCK TO PREPARE LOG LINES.
C
C///////////////////////////////////////////////////////////////////////

7100  CONTINUE

      DO 7200 J = 1, 6
         COLUMN(J) = ' '
7200  CONTINUE

      REMARK = ' '

C     COLUMN 1: NAME OF THE ACTIVITY
      IF (QTASK .EQ. QENTRY) COLUMN(1) = '    START'
      IF (QTASK .EQ. QNEWTO) COLUMN(1) = '   NEWTON'
      IF (QTASK .EQ. QREFIN) COLUMN(1) = '   REFINE'
      IF (QTASK .EQ. QTIMST) COLUMN(1) = '   TIMSTP'

C     COLUMN 2: NORM OF THE STEADY-STATE FUNCTION
      IF (.NOT. FOUND) GO TO 7400
         ASSIGN 7300 TO RETURN
         GO TO 9941
7300     CONTINUE
         CALL NORM (SCALRS + COMPS * POINTS, TEMP, BUFFER)
         CALL LOGSTR ('(F9.2)', 9, COLUMN(2), TEMP)
7400  CONTINUE

C     COLUMN 3: LARGEST CONDITION NUMBER
      IF (QTASK .EQ. QNEWTO .OR. QTASK .EQ. QTIMST) THEN
         IF (MAXCON .NE. ZERO)
     +      CALL LOGSTR ('(F9.2)', 9, COLUMN(3), MAXCON)
      END IF

C     COLUMN 4: NUMBER OF POINTS
      IF (QTASK .EQ. QENTRY .OR. QTASK .EQ. QREFIN) THEN
         WRITE (COLUMN(4), '(I4)') POINTS
      END IF

C     COLUMN 5: NUMBER OF STEPS
      IF (QTASK .EQ. QNEWTO .OR. QTASK .EQ. QTIMST) THEN
         WRITE (COLUMN(5), '(I4)') STEPS
      END IF

C     COLUMN 6: NUMBER OF JACOBIANS
      IF (QTASK .EQ. QNEWTO .OR. QTASK .EQ. QTIMST) THEN
         WRITE (COLUMN(6), '(I4)') JACOBS
      END IF

C     REMARK
      IF (QTASK .EQ. QNEWTO .OR. QTASK .EQ. QTIMST) THEN
         IF (XREPOR .EQ. QDVRG) THEN
            REMARK = 'DIVERGING'
         ELSE IF (XREPOR .EQ. QNULL) THEN
            REMARK = ' '
         ELSE IF (XREPOR .EQ. QBNDS) THEN
            REMARK = 'AT BOUNDS'
         ELSE IF (XREPOR .EQ. QDONE) THEN
            REMARK = ' '
         ELSE IF (XREPOR .EQ. QSMAX) THEN
            REMARK = 'STEP LIMIT'
         ELSE IF (XREPOR .EQ. QSTDY) THEN
            REMARK = 'STEADY STATE'
         ELSE
            REMARK = '?'
         END IF
      ELSE IF (QTASK .EQ. QREFIN) THEN
         IF (FOUND) REMARK = 'NEW MESH'
      END IF

      GO TO LABEL

C///////////////////////////////////////////////////////////////////////
C
C     REQUEST REVERSE COMMUNICATION.
C
C///////////////////////////////////////////////////////////////////////

C///  SAVE THE SOLUTION.

9911  CONTINUE
      CALL CPUTIM (TIMER(QOTHER))

      CALL COPY (SCALRS + COMPS * POINTS, X, BUFFER)
      SAVE = .TRUE.
      REENT = .TRUE.
      ASSIGN 9912 TO ROUTE
      GO TO 99999
9912  CONTINUE

      CALL ELAPSE (TIMER(QOTHER))
      TOTAL(QOTHER) = TOTAL(QOTHER) + TIMER(QOTHER)
      IF (GRID .LE. GMAX) THEN
         DETAIL(GRID, QOTHER) = DETAIL(GRID, QOTHER) + TIMER(QOTHER)
         EVENT(GRID, QOTHER) = EVENT(GRID, QOTHER) + 1
      END IF
      GO TO RETURN

C///  PRINT THE LATEST SOLUTION.

9921  CONTINUE
      CALL CPUTIM (TIMER(QOTHER))

      CALL COPY (SCALRS + COMPS * POINTS, X, BUFFER)
      SHOW = .TRUE.
      REENT = .TRUE.
      ASSIGN 9922 TO ROUTE
      GO TO 99999
9922  CONTINUE

      CALL ELAPSE (TIMER(QOTHER))
      TOTAL(QOTHER) = TOTAL(QOTHER) + TIMER(QOTHER)
      IF (GRID .LE. GMAX) THEN
         DETAIL(GRID, QOTHER) = DETAIL(GRID, QOTHER) + TIMER(QOTHER)
         EVENT(GRID, QOTHER) = EVENT(GRID, QOTHER) + 1
      END IF
      GO TO RETURN

C///  PASS REQUESTS FROM NEWTON, REFINE, OR TIMSTP TO THE CALLER.

9931  CONTINUE

C     IDENTIFY THE REQUEST.  THIS MUST BE SAVED TO GATHER STATISTICS
C     AT REENTRY.  THE REVERSE COMMUNICATION FLAGS WILL NOT BE SAVED
C     BECAUSE THEY ARE CLEARED AT EVERY ENTRY.
      IF (FUNCT) THEN
         QTYPE = QFUNCT
      ELSE IF (JACOB) THEN
         QTYPE = QJACOB
      ELSE IF (SOLVE) THEN
         QTYPE = QSOLVE
      ELSE
         QTYPE = QOTHER
      END IF

C     COUNT THE NUMBER OF JACOBIANS FOR INCLUSION IN THE LOG LINE.
      IF (JACOB) JACOBS = JACOBS + 1

      CALL CPUTIM (TIMER(QTYPE))

      REENT = .TRUE.
      ASSIGN 9932 TO ROUTE
      GO TO 99999
9932  CONTINUE

C     SAVE THE MAXIMUM CONDITION NUMBER.
      IF (QTYPE .EQ. QJACOB) MAXCON = MAX (MAXCON, CONDIT)

      CALL ELAPSE (TIMER(QTYPE))
      TOTAL(QTYPE) = TOTAL(QTYPE) + TIMER(QTYPE)
      IF (GRID .LE. GMAX) THEN
         DETAIL(GRID, QTYPE) = DETAIL(GRID, QTYPE) + TIMER(QTYPE)
         EVENT(GRID, QTYPE) = EVENT(GRID, QTYPE) + 1
      END IF
      GO TO RETURN

C///  EVALUATE THE STEADY-STATE FUNCTION.

9941  CONTINUE
      CALL CPUTIM (TIMER(QFUNCT))

      CALL COPY (SCALRS + COMPS * POINTS, X, BUFFER)
      FUNCT = .TRUE.
      TIME = .FALSE.
      REENT = .TRUE.
      ASSIGN 9942 TO ROUTE
      GO TO 99999
9942  CONTINUE

      CALL ELAPSE (TIMER(QFUNCT))
      TOTAL(QFUNCT) = TOTAL(QFUNCT) + TIMER(QFUNCT)
      IF (GRID .LE. GMAX) THEN
         DETAIL(GRID, QFUNCT) = DETAIL(GRID, QFUNCT) + TIMER(QFUNCT)
         EVENT(GRID, QFUNCT) = EVENT(GRID, QFUNCT) + 1
      END IF
      GO TO RETURN

C///////////////////////////////////////////////////////////////////////
C
C     INFORMATIVE MESSAGES.
C
C///////////////////////////////////////////////////////////////////////

10001 FORMAT
     +  (/1X, A9, A
     +   /10X, A
     +  //20X, A
     +   /20X, A
     +   /20X, A)

10002 FORMAT
     +  (/1X, A9, 'INITIAL GUESS:')

10003 FORMAT
     +  (/1X, A9, 'LOG OF THE BOUNDARY VALUE PROBLEM SOLVER.'
     +  //10X, A36, A26
     +   /10X, A36, A26
     +   /10X, A36, A26)

10004 FORMAT
     +  (/10X, A9,
     +    3X, A9,
     +    3X, A9,
     +    3X, A4,
     +    1X, A4,
     +    1X, A4,
     +    3X, A9,
     +    3X, A9)

10005 FORMAT
     +  (/1X, A9, 'FINAL SOLUTION:')

10006 FORMAT
     +  (/1X, A9, 'TOTAL CPU TIME.'
     +  //10X, A9, 2X, A)

10007 FORMAT
     +  (/1X, A9, 'PERCENT OF TOTAL CPU TIME.'
     +   /3(/10X, A16, A22, A27)
     +  //(10X, I6, 2X, F6.1, 1X, 3(1X, F6.1), 1X, 4(1X, F6.1)))

10008 FORMAT
     +  (/10X, ' TOTAL', 9X, 3(1X, F6.1), 1X, 4(1X, F6.1))

10009 FORMAT
     +  (/10X, 'SPACE LIMITATIONS PREVENT RETAINING TIMING DATA FOR'
     +   /10X, 'EVERY GRID.  NEVERTHELESS, THE TOTALS REFER TO ALL'
     +   /10X, 'GRIDS; THEY ARE NOT COLUMN SUMS.')

10010 FORMAT
     +  (/1X, A9, 'AVERAGE CPU TIME.'
     +   /3(/10X, A9, A28, A25)
     +  //(10X, I6, 3X, F7.3, 2X, F7.2, 2X, F7.3, 1X, 3(2X, I7)))

10011 FORMAT
     +  (/10X, 'SPACE LIMITATIONS PREVENT RETAINING TIMING DATA FOR'
     +   /10X, 'EVERY GRID.')

10012 FORMAT
     +  (/1X, A9, 'SUCCESS.  BOUNDARY VALUE PROBLEM SOLVED.')

10013 FORMAT
     +  (/1X, A9, 'FAILURE.  NO SOLUTION WAS FOUND, NOT EVEN ON THE'
     +   /10X, 'FIRST GRID.')

10014 FORMAT
     +  (/1X, A9, 'FAILURE.  A SOLUTION WAS FOUND FOR SOME GRIDS, BUT'
     +   /10X, 'NOT THE LAST.')

10015 FORMAT
     +  (/1X, A9, 'FAILURE.  THE SOLUTION IS NOT RESOLVED TO THE EXTENT'
     +   /10X, 'DESIRED BECAUSE THERE IS NO SPACE FOR MORE POINTS.')

10016 FORMAT
     +  (/1X, A9, 'CALLING NEWTON TO SOLVE THE STEADY-STATE PROBLEM.')

10017 FORMAT
     +  (10X, A9,
     +   3X, A9,
     +   3X, A9,
     +   3X, A4,
     +   1X, A4,
     +   1X, A4,
     +   3X, A)

10018 FORMAT
     +  (/1X, A9, 'NEWTON SOLVED THE STEADY-STATE PROBLEM.')

10019 FORMAT
     +  (/1X, A9, 'NEWTON DID NOT SOLVE THE STEADY-STATE PROBLEM.')

10020 FORMAT
     +  (/1X, A9, 'CALLING REFINE TO ADAPT THE MESH.')

10021 FORMAT
     +  (/1X, A9, 'REFINE PRODUCED A NEW MESH.')

10022 FORMAT
     +  (/1X, A9, 'REFINE DID NOT PRODUCE A NEW MESH.')

10023 FORMAT
     +  (/1X, A9, 'CALLING TIMSTP TO PERFORM A TIME INTEGRATION.')

10024 FORMAT
     +  (/1X, A9, 'TIMSTP COMPLETED THE TIME INTEGRATION.')

10025 FORMAT
     +  (/1X, A9, 'TIMSTP DID NOT COMPLETE THE TIME INTEGRATION.')

C///////////////////////////////////////////////////////////////////////
C
C     ERROR MESSAGES.
C
C///////////////////////////////////////////////////////////////////////

C     GO TO 99999

9001  IF (0 .LT. TEXT) WRITE (TEXT, 99001) ID
      GO TO 99999

9002  IF (0 .LT. TEXT) WRITE (TEXT, 99002) ID, PASSES
      GO TO 99999

9003  IF (0 .LT. TEXT) WRITE (TEXT, 99003) ID,
     +   SCALRS, COMPS, POINTS, PMAX
      GO TO 99999

9004  IF (0 .LT. TEXT) WRITE (TEXT, 99004) ID,
     +   COMPS, PMAX, SCALRS + COMPS * PMAX
      GO TO 99999

9005  IF (0 .LT. TEXT) WRITE (TEXT, 99005) ID
      GO TO 99999

9006  IF (0 .LT. TEXT) WRITE (TEXT, 99006) ID, LAST, ISIZE
      GO TO 99999

9007  IF (0 .LT. TEXT) WRITE (TEXT, 99007) ID, LAST, RSIZE
      GO TO 99999

9008  IF (0 .LT. TEXT) WRITE (TEXT, 99008) ID
      GO TO 99999

9009  IF (0 .LT. TEXT) WRITE (TEXT, 99009) ID
      GO TO 99999

9010  IF (0 .LT. TEXT) WRITE (TEXT, 99010) ID
      GO TO 99999

9011  IF (0 .LT. TEXT) WRITE (TEXT, 99011) ID
      GO TO 99999

9012  IF (0 .LT. TEXT) WRITE (TEXT, 99012) ID
      GO TO 99999

99001 FORMAT
     +  (/1X, A9, 'ERROR.  THE PRINTING LEVELS ARE NOT ORDERED.'
     +  //10X, I10, '  LEVEL(1) FOR MESSAGES'
     +   /10X, I10, '  LEVEL(2) FOR THE SOLUTION')

99002 FORMAT
     +  (/1X, A9, 'ERROR.  THE NUMBER OF PASSES MUST BE POSITIVE.'
     +  //10X, I10, '  PASSES')

99003 FORMAT
     +  (/1X, A9, 'ERROR.  THE NUMBERS OF VARIABLES ARE INCONSISTENT.'
     +  //10X, I10, '  SCALARS'
     +   /10X, I10, '  COMPONENTS'
     +   /10X, I10, '  POINTS'
     +   /10X, I10, '  LIMITING NUMBER OF POINTS')

99004 FORMAT
     +  (/1X, A9, 'ERROR.  THE ARRAY DIMENSIONS ARE MUST BE POSITIVE.'
     +  //10X, I10, '  COMPS'
     +   /10X, I10, '  PMAX'
     +   /10X, I10, '  SCALRS + COMPS * PMAX')

99005 FORMAT
     +  (/1X, A9, 'ERROR.  GRAB FAILS.')

99006 FORMAT
     +  (/1X, A9, 'ERROR.  THE INTEGER WORK SPACE IS TOO SMALL.'
     +   //10X, I10, '  NEEDED'
     +    /10X, I10, '  AVAILABLE')

99007 FORMAT
     +  (/1X, A9, 'ERROR.  THE REAL WORK SPACE IS TOO SMALL.'
     +   //10X, I10, '  NEEDED'
     +    /10X, I10, '  AVAILABLE')

99008 FORMAT
     +  (/1X, A9, 'ERROR.  UNKNOWN TASK.')

99009 FORMAT
     +  (/1X, A9, 'ERROR.  UNKNOWN REPORT CODE.')

99010 FORMAT
     +  (/1X, A9, 'ERROR.  NEWTON ERRS.')

99011 FORMAT
     +  (/1X, A9, 'ERROR.  REFINE ERRS.')

99012 FORMAT
     +  (/1X, A9, 'ERROR.  TIMSTP ERRS.')

C///  EXIT.

99999 CONTINUE
      RETURN
      END
